矩阵数值分析上机作业

《矩阵与数值分析》课程数值实验报告

教学班号：

任课老师：

学部：

专业：

姓名：

学号：

1.方程在x=3.0附近有根，试写出其三种不同的等价形

式以构成两种不同的迭代格式，再用这两种迭代求根，并绘制误差下降曲线，观察这两种迭代是否收敛及收敛的快慢

newton迭代法

function x=juz1(x1,tol,max)

%x1为初值，tol为误差终止条件，max为最大迭代次数

xo=x1;

T=[];%创建T矩阵储存每次迭代的结果

i=1;

while i

xn=xo-(2*xo^3-5*xo^2-19*xo+42)/(6*xo^2-10*xo-19);

T(i)=xn;

if abs(xn-xo)

disp('迭代法收敛');

disp([xn]);

R=xn-T;%计算误差

a=linspace(0,0.1,i);

plot(a,R,'r');%绘制误差下降曲线

title('误差下降曲线');

return;

else

xo=xn;

i=i+1;

end

disp('迭代法不收敛');

disp('迭代最大次数的结果为：');

>> juz1(3,10^-5,20)

迭代法收敛

3.5000

i =

弦截法

function x=juz11(x1,x2,tol,max)

%x1,x2为初值，tol为误差终止条件，max为最大迭代次数

xo=x1;

xp=x2;

i=1;

while i

xn=xp-(2*xp^3-5*xp^2-19*xp+42)/(2*xp^2+2*xp*xo+2*xo^2-5*xp-5*xo-19); T(i)=xn;

if abs(xn-xp)

disp('迭代法收敛');

disp([xn]);

R=xn-T;

a=linspace(0,1,i);

plot(a,R,'r');%绘制误差下降曲线

title('误差下降曲线');

return;

else

xo=xp;

xp=xn;

i=i+1;

end

disp('迭代法不收敛');

disp('迭代最大次数的结果为：');

>> juz11(2.9,3,10^-5,20)

迭代法收敛

3.5000

i =

2.用复化梯形公式、复化辛普森公式、龙贝格公式求下列定积分，要求绝对误差为

，并将计算结果与精确解进行比较：

（1）（2）

（1）用复化梯形公式求积分

function f=Txing1(a,b,n)

% a，b分别为积分下限和上限，n为划分的区间个数

while n<10^4

h1=(b-a)/n;%h1为积分区间长度

x1=zeros(1,n+1);

for i=1:n+1

x1(i)=a+(i-1)*h1;%计算n+1个节点储存在x1中

end

y1=(2/3)*(x1.^3).*exp(x1.^2);%计算每个节点对应的函数值 to=0;

for i=1:n

to=to+h1/2*(y1(i)+y1(i+1));%用梯形公式求积分

end

h2=(b-a)/(n+1);%h2为积分区间长度

x2=zeros(1,n+2);

for i=1:n+2

x2(i)=a+(i-1)*h2;%计算n+2个节点储存在x2中

end

y2=(2/3)*(x2.^3).*exp(x2.^2);%计算每个节点对应的函数值 tn=0;

for i=1:n+1

tn=tn+h2/2*(y2(i)+y2(i+1));%用梯形公式求积分

end

if abs(tn-to)<0.5*10^(-8)%进行误差判断

I=tn %计算结果

N=n+1 %实际分割区间个数

r=exp(4)-tn %计算误差

return;

else

n=n*2;

end

disp('分割区间太长达不到精度要求');

n/2+1

>> Txing1(1,2,5)

I =

54.5982

N =

5121

r =

-5.0604e-06

复化辛普森求积分

function f=Simpson1(a,b,n)

%a,b分别为积分下限和上限

%n为分割区间个数

while n<10^8

T=[];

for n=n:n+1

h1=(b-a)/n;%将积分区间n等分

x1=zeros(1,n+1);

for i=1:n+1

x1(i)=a+(i-1)*h1;%计算n+1个求积节点值并储存在x中

end

z1=(2/3)*x1.^3.*exp(x1.^2);%计算n+1个求积节点处的函数值

y1=zeros(1,n);

for i=1:n

y1(i)=a+(2*i-1)*h1/2;%计算n个区间中点值

end

r1=(2/3)*y1.^3.*exp(y1.^2);%计算n个区间中点对应的函数值

t=0;

for i=1:n;

t=t+h1/6*(z1(i)+4*r1(i)+z1(i+1));%用复化simpson公式求积分end

T=[T,t];

end

if abs(T(2)-T(1))<0.5*10^(-8)%进行误差判断

I=T(2) %计算结果

N=n+1 %实际分割区间个数

r=exp(4)-T(2) %计算误差

return;

else

n=2*n;

end

disp('分割区间太长达不到精度要求');

T(2)

n/2+1

>> Simpson1(1,2,3)

I =

54.5982

N =

160

r =

-2.9586e-08

龙贝格公式求积分

function f=romberg1(a,b)

t(1,1)=(b-a)/2*(2*a^3*exp(a^2)/3+2*b^3*exp(b^2)/3);

n=1;

while n<100

for k=1:n

aa=0;

for i=0:(power(2,(k-1))-1);

aa=aa+2*(a+(2*i+1)*(b-a)/2^k)^3*exp((a+(2*i+1)*(b-a)/2^k)^2)/3; end

t(1,k+1)=0.5*t(1,k)+(b-a)/(2^k)*aa;

end

for m=1:n

h=1/((4^m)-1);

for k=m:n

t(m+1,k+1)=((4^m)*t(m,k+1)-t(m,k))*h;

end

if (abs(t(n,n)-t(n+1,1+n))<0.5*10^(-8))

I=t(n,n)

return;

else

n=n+1;

end

disp('超出最大迭代次数');

>> romberg1(1,2)

I =

54.5982

（2）复化梯形公式求积分

function f=Txing2(a,b,n)

% a，b分别为积分下限和上限，n为划分的区间个数

while n<10000

h1=(b-a)/n;%h为积分区间长度

x1=zeros(1,n+1);

for i=1:n+1

x1(i)=a+(i-1)*h1;%计算n+1个节点储存在x1中 end

y1=2*x1./(x1.^2-3);%计算每个节点对应的函数值

to=0;

for i=1:n

to=to+h1/2*(y1(i)+y1(i+1));%用梯形公式求积分 end

h2=(b-a)/(n+1);%h2为积分区间长度

x2=zeros(1,n+2);

for i=1:n+2

x2(i)=a+(i-1)*h2;%计算n+2个节点储存在x2中 end

y2=2*x2./(x2.^2-3);%计算每个节点对应的函数值

tn=0;

for i=1:n+1

tn=tn+h2/2*(y2(i)+y2(i+1));%用梯形公式求积分 end

if abs(tn-to)<0.5*10^(-8)%进行误差判断

I=tn %计算结果

N=n+1 %实际分割区间个数

r=log(6)-tn %计算结果与精确值比较 return;

else

n=n*2;

end

disp('分割区间太长达不到精度要求');

n/2+1

>> Txing2(2,3,4)

I =

1.7918

N =

1025

r =

-1.0576e-06

复化辛普森公式求积分

function f=Simpson2(a,b,n)

%a,b分别为积分下限和上限

%n为分割区间个数

while n<10000

T=[];%生成空矩阵以储存计算结果

for n=n:n+1

h=(b-a)/(2*n);%将积分区间分成2n份

x=zeros(1,2*n+1);

for i=1:2*n+1

x(i)=a+(i-1)*h;%计算n+1个求积节点值并储存在x中

end

z=2*x./(x.^2-3);%计算2n+1个求积节点处的函数值

t=0;

for i=1:n;

t=t+h/3*(z(2*i-1)+4*z(2*i)+z(2*i+1));%用复化simpson公式求积分end

T=[T,t];%将分成n份和n+1份结果储存在T中

end

if abs(T(2)-T(1))<0.5*10^(-8)%进行误差判断

I=T(2) %计算结果

N=n+1 %实际分割区间个数

r=log(6)-t %计算误差

return;

else

n=2*n;

end

disp('分割区间太长达不到精度要求');

T(2)

n/2+1

>> Simpson2(2,3,4)

I =

1.7918

N =

r =

-4.9476e-09

龙贝格公示求积分

function f=romberg2(a,b)

t(1,1)=(b-a)/2*(2*a/(a^2-3)+2*b/(b^2-3));

n=1;

while n<100

for k=1:n

aa=0;

for i=0:(power(2,(k-1))-1);

aa=aa+2*(a+(2*i+1)*(b-a)/2^k)/((a+(2*i+1)*(b-a)/2^k)^2-3); end

t(1,k+1)=0.5*t(1,k)+(b-a)/(2^k)*aa;

end

for m=1:n

h=1/((4^m)-1);

for k=m:n

t(m+1,k+1)=((4^m)*t(m,k+1)-t(m,k))*h;

end

if (abs(t(n,n)-t(n+1,1+n))<0.5*10^(-8))

I=t(n,n)

return;

else

n=n+1;

end

disp('超出最大迭代次数');

>> romberg2(2,3)

I =

1.7918

3.使用带选主元的分解法求解线性方程组,其中，，

当时．对于的情况分别求解．

精确解为．对得到的结果与精确解的差异进行解释．

function f=LU(n)

%生成矩阵A

A=zeros(n);

for i=1:n

for j=1:n

A(i,j)=i^(j-1);

end

%生成右边项b

b=zeros(n,1);

b(1)=n;

for i=2:n

b(i)=(i^n-1)/(i-1);

end

x=zeros(n,1);

y=zeros(n,1);

a=zeros(n,1);%储存矩阵A的换行元素

c=0;%储存右边项的换行元素

for i=1:n-1

[e,m]=max(abs(A(i:n,i)));%寻找最大元素及其位置

a=A(i,:);

A(i,:)=A(m+i-1,:);

A(m+i-1,:)=a;

c=b(i);

b(i)=b(m+i-1);

b(m+i-1)=c;

if A(i,i)==0

disp('A是奇异矩阵，方程组没有唯一解');

return;

end

for j=i+1:n

d=A(j,i)/A(i,i);

A(j,i)=d;

A(j,i+1:n)=A(j,i+1:n)-d*A(i,i+1:n);

end

%求解Ly=b

y(1)=b(1);

for k=2:n

y(k)=b(k)-A(k,1:k-1)*y(1:k-1);

end

%求解Ux=y x(n)=y(n)/A(n,n); for k=n-1:-1:1

x(k)=(y(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end x

>> LU(3) x =

1 1 1 >> LU(7) x =

1 1 1 1 1 1 1 >> LU(11) x =

1.0001 0.9997 1.0004 0.9998 1.0001 1.0000 1.0000 1.0000 1.0000 1.0000

1.0000

进行LU 分解时，系数矩阵出现病态，产生误差 4. 用4阶Runge-kutta 法求解微分方程

t t t te e t u u u u u 22210

)(,101)0(,2---+==

-=' (1) 令1.0=h ，使用上述程序执行20步，然后令05.0=h ，使用上述程序执行40步

(2) 比较两个近似解与精确解

(3) 当h减半时，(1)中的最终全局误差是否和预期相符？

(4) 在同一坐标系上画出两个近似解与精确解．（提示输出矩阵R包含近似解的x和y坐

标，用命令plot(R(:,1),R(:,2))画出相应图形．）

function f=Rungek(x,t)

uo=x;

tn=t;

%h=0.1,计算20步，并将计算结果储存在U1中

U1=zeros(20,1);

h1=0.1;

i=1;

while i<=20

k1=exp(-2*tn)-2*uo;

k2=exp(-2*(tn+0.5*h1))-2*(uo+0.5*h1*k1);

k3=exp(-2*(tn+0.5*h1))-2*(uo+0.5*h1*k2);

k4=exp(-2*(tn+h1))-2*(uo+h1*k3);

un=uo+h1/6*(k1+2*k2+2*k3+k4);

U1(i)=un;

i=i+1;

tn=(i-1)*h1;

uo=un;

end

%h=0.05,计算40步，并将计算结果储存在U2中

U3=zeros(40,1);

h2=0.05;

i=1;

uo=x;

tn=t;

while i<=40

s1=exp(-2*tn)-2*uo;

s2=exp(-2*(tn+0.5*h2))-2*(uo+0.5*h2*s1);

s3=exp(-2*(tn+0.5*h2))-2*(uo+0.5*h2*s2);

s4=exp(-2*(tn+h2))-2*(uo+h2*s3);

un=uo+h2/6*(s1+2*s2+2*s3+s4);

U3(i)=un;

i=i+1;

tn=(i-1)*h2;

uo=un;

end

U2=zeros(20,1);

for i=1:20

U2(i)=U3(2*i);

end

%利用函数U的表达式计算出精确结果，并储存在U中

U=zeros(20,1);

for i=1:20

t=i*0.1;

U(i)=0.1*exp(-2*t)+t*exp(-2*t);

end

%显示计算结果，其中R1是U1的绝对误差，R2是U2的绝对误差，R3是U1与U2的差U1,U2,U,R1=U-U1,R2=U-U2,R3=U1-U2

a=zeros(20,1);

for i=1:20

a(i)=i;

end

plot(a,U1,'+',a,U2,'o',a,U2,'r');

>> Rungek(0.1,0)

U1 =

0.1637

0.2011

0.2195

0.2247

0.2207

0.2108

0.1973

0.1817

0.1653

0.1489

0.1330

0.1179

0.1040

0.0912

0.0797

0.0693

0.0601

0.0519

0.0447

0.0385

U2 =

0.1637

0.2011

0.2195

0.2247

0.2207

0.2108

0.1973

0.1817

0.1653

0.1489

0.1330

0.1179

0.1040

0.0912

0.0797

0.0693

0.0601

0.0519

0.0447

0.0385 U =

0.1637

0.2011

0.2195

0.2247

0.2207

0.2108

0.1973

0.1817

0.1653

0.1489

0.1330

0.1179

0.1040

0.0912

0.0797

0.0693

0.0601

0.0519

0.0447

0.0385 R1 =

1.0e-005 *

0.2114

0.3250

0.3732

0.3791

0.3590

0.3242

0.2825

0.2389

0.1966

0.1575

0.1226

0.0924

0.0667

0.0454

0.0281

0.0142

0.0034

-0.0048

-0.0108

-0.0151 R2 =

1.0e-006 *

0.1191

0.1830

0.2098

0.2127

0.2010

0.1811

0.1574

0.1326

0.1087

0.0866

0.0670

0.0499

0.0356

0.0236

0.0140

0.0063

0.0003

-0.0042

-0.0075

-0.0097 R3 =

1.0e-005 *

-0.1995

-0.3067

-0.3522

-0.3578

-0.3389

-0.3061

-0.2667

-0.2256

-0.1857

-0.1488

-0.1159

-0.0874

-0.0632

-0.0430

-0.0267

-0.0136

-0.0034

0.0044

0.0101

0.0141

5.设

为阶的三对角方阵，是一个阶的对称正定矩阵

其中为阶单位矩阵。设为线性方程组的真解，右边的向量由这个真解给出。

(1)用Cholesky分解法求解该方程.

(2)用Jacobi迭代法和Gauss-Seidel迭代法求解该方程组，误差设为.

其中取值为4，5，6.

用Cholesky分解法求解方程：

function f=Cholesky(n)

%生成Tn矩阵

T=zeros(n);

T(i,i)=2;

end

for i=1:n-1

T(i,i+1)=-1;

end

for i=2:n

T(i,i-1)=-1;

end

%生成系数矩阵A

I=eye(n);

B=T+2*I;

C=-I;

A=zeros(n^2);

for i=1:n

A(n*(i-1)+1:i*n,n*(i-1)+1:i*n)=B; end

for i=1:n-1

A(i*n+1:(i+1)*n,(i-1)*n+1:i*n)=C; A((i-1)*n+1:i*n,i*n+1:(i+1)*n)=C; end

%生成右边项b

y=ones(n^2,1);

b=A*y;

R=chol(A);%对系数矩阵A进行Cholesky分解x=R\(R'\b)%求解方程组

>> Cholesky(4)

x =

1.0000

>> Cholesky(5) x =

1.0000

>> Cholesky(6) x =

1.0000