2023年矩阵实验报告

矩阵与数值分析实验报告

任课教师:________张宏伟_________________________

教学班号：__________02___________________________

院系：电信（计算机应用技术）

学生姓名:

学生学号:

1.方程,f(c)=2ir3-5rr2-19a:+42=0在x=3.0附近有根,试写出其三种不同的

等价形式以构成两种不同的迭代格式，再用这两种迭代求根，并绘制误差下

降曲线，观测这两种迭代是否收敛及收敛的快慢。

解答：

代码如下：

clear;

symsxty；

x=3;%初始迭代点

t=0;%中间变量

y=o；%绘制下降曲线变化量

k=0；/迭代计数变量

epx=1;%变量计算差

E=le-20;%精度

f1=(2*x~3-5*x-2+42)/19;%迭代1

f2=(2*x-3-19*x+42)”(l/2);%迭代2

f3=(5*x*2+19*x-42)*(1/3)；%迭代3

whi1e(epx>E&k<1000)%循环迭代

k=k+l;

y(k)=x；

t=f1；%标记1

epx=abs(t-x);

x=t;

f1=(2*x*3-5*x,2+42)/19；%标记2

end；

Plot(y);

title（'迭代1误差下降曲线'）;

迭代公式1收敛结果：x=2

迭代公式1误差变化曲线

x10迭代2误差下降曲线

3.5

2.5

2

1.5

0.5

10

迭代公式2误差变化曲线

迭代公式3收敛结果:x=6.8750

迭代公式3误差变化曲线

结果分析:迭代公式1：fl=(2*x*3-5*x*2+42)/19和迭代公式3：f3=(5*x-2

+19*x-42)71/3)计算结果收敛,其中迭代公式1收敛速度快于迭代公式3o

迭代公式2：f2=(2*x*3-19*x+42)71/2)计算结果不收敛。

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

差为,=().5x10-8并将计算结果与精确解进行比较：

4232

e=-£e.2n6=

3

解答：

(1)复化梯形公式代码：

clear;

symsxsum1sum2albl;

fun=(2/3)*(x-3)*exp(x-2);

suml=0;%积分结果变量1

sum2=0；%积分结果变量2

n=1；%迭代计数变量

epx=l阶段误差

E=0.5e-8;%精度

a=l;%积分下限

b=2;%积分上限

while(epx>E)

h=(b—a)/n；

fori=0:(nT)%for循环求解一次N点积分结果

a1=subs(fun,x,(a+i*h))；

bl=subs(fun,x,(a+(i+1)*h))；

t=((al+bl)*h)/2；

sum2=sum2+t；

end;

epx=abs(suml-sum2);%计算阶段误差

suml=sum2;为使用suml暂存上次的计算结果

sum2=0；

n=n*2;

end；

disp('积分结果为：，),

vpa(sum1)

积分结果为：

ans=

54.4316975

真值约为：54.598

分析：由结果可知计算结果具有8位有效数字，已满足精度规定。使用数据点数

为N=8192o

复化辛普森公式代码：

clear;

symsxsumlsum2a1b1;

fun=(2/3)*(x-3)*exp(x'2);

suml=0;%积分结果变量1

sum2=0；%积分结果变量2

n=l;%迭代计数变量

epx=1;%变量计算差

E=0.5e-8；%精度

a=l;%积分下限

b=2;%积分上限

whi1e(epx>E)

h=(b—a)/n;

fori=0:(n-1)

a1=subs(fun,x,(a+i*h));

bl=subs(fun,x,(a+(i+1)*h))；

ab=subs(fun,x,((a+i*h)+(a+(i+l)*h))/2);

t=((al+4*ab+bl)*h)/6；

sum2=sum2+t；

end；

epx=abs(sum1-sum2)；

suml=sum2;

sum2=0;

n=n*2;

end；

disp('积分结果为:’)，

vpa(suml)

积分结果为：

ans=

54.59858433

真值约为：54.598

分析：由结果可知计算结果具有11位有效数字，已满足精度规定。使用数据点

数为N=1024,对比复化梯形公式可以复化辛普森公式计算精度更高。

龙贝格公式代码：

clear；

symsxa1bly;

fun=(2/3)*(x-3)*exp(x-2);

suml=O;96积分结果变量1

sum2=0;%积分结果变量2

epx=l;%变量计算差

E=0.5e-8；蟒青度

a=1;%积分下限

b=2;%积分上限

k=0;%迭代计数变量1

m=0;%迭代计数变量2

while(epx>E)

k=k+1;

m=m+1;

h=(b—a)/(2-k);

fori=0:⑵k—l)

al=subs(fun,x,(a+i*h));

b1=subs(fun,x,(a+(i+1)*h));

t=((al+bl)*h)/2;

sum2=sum2+t；

end;

sum1=sum2；

y(m)=sum2;

sum2=0;

fori=0:k-2;

m=m+1；

y(m)=((4"(k—1))*y(m—1)—y(m—k))/(4^(k-l)—1);

end；

if(m>l)

epx=abs(sym(y(m)-y(m-1)));

end;

end;

disp(,积分结果为：'),

Vpa(y(m))

积分结果为：

ans=

54.6804633

真值约为：54.598

分析：由结果可知计算结果具有9位有效数字，已满足精度规定。

(2)积分结果为：

ans=

1.96865767

真值约为：1.

分析：由结果可知计算结果具有8位有效数字，已满足精度规定。使用数据点数为

N=65536O

积分结果为：

ans=

1.48476064

真值约为：1.

分析：由结果可知计算结果具有10位有效数字，已满足精度规定。使用数据点数

为N=512,可知复化辛普森公式精度明显高于复化梯形公式。

积分结果为：

ans=

1.

真值约为：1.

分析：由结果可知计算结果具有8位有效数字，已满足精度规定。

3.使用带选主元的分解法求解线性方程组AT=h,其中4=S’lnx"，3=",

^ilnxi,bn=n.当，＞2时,bn=(*-1)/("1).对于"=37”的情况分别求解。

精确解为工=(11---IIf对得到的结果与精确解的差异进行解释。

解答：

程序代码：

functionx=luso1ve(A,b)

[n,n]=size(A);

x=zeros(n,1);

y=zeros(n,1);

temprow=zeros(n,1)；

tempconstant=0;

Pvector=zeros(n,1);

forcol=l:n-1

[max_element,index]=max(abs(A(co1:n,col)));

temprow=A(col,:)；

A(col,：)=A(index+co1-1,:);

A(index+col-1,:)=temprow;

tempconstant=b(co1);

b(col)二b(index+co1—1);

b(index+col-1)=tempconstant;

ifA(co1,col)=0

disp('Aissingular.nouniquesolution');

return

end

forrow=col+1:n

mult=A(row,col)/A(col,col);

A(row,col)=mult;

A(row,col+1:n)=A(row,col+1:n)—mul

t*A(col,col+1:n)；

end

end

y(1)=b(l)；

fork=2：n

y(k)=b(k)—A(k,1：k-1)*y(l：k-1);

end

x(n)=y(n)/A(n,n);

fork=n-1：-l:1

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

end

实验结果：

T

n=3时,解得x=[1,1,1];

T

n=7时，解得x=[1,1,1,1,1>1,1]；

n=11时，解得x=[1.0001,0.9998,1.0002,0.9999,1,1,1,1,1,

1,I；

实验结果分析：

高斯列主元消去法，避免了小主元做除数，在Gauss消去法中增长选主

元的过程,在第K步消元时,一方面在第K列主对角元以下元素中挑选绝对值

最大的数，并通过初等行变换，使得该数位于主对角上，然后继续消元。当矩

阵维数较小时，精度较高，但随着矩阵维数增大,精度下降。

4.用4阶Runge-kutta法求解微分方程

u=e~2'-2u,u(0)=—,u(t)=-e-%+te2'

1010

(1)令也=0］，使用上述程序执行20步，然后令人=005,使用上述程序执行4

0步

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

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

(4)在同一坐标系上画出两个近似解与精确解.(提醒输出矩阵R包含近似解的

x和V坐标，用命令pl。t(R(:,1),R(:,2))画出相应图形.)

解答：

程序代码：

functionR=rk4(f,a,b,ya,n)

symstu

h=(b—a)/n;

T=zeros(1,n+1);

Y=zeros(l,n+1)；

T=a:h:b;

Y(1)=ya;

fork=1:n

K1=h*subs(f,{t,u},{T(k),Y(k)});

K2=h*subs(f,{t,u},{T(k)+h/2,Y(k)+Kl/2});

K3=h*subs(f,{t,u},{T(k)+h/2,Y(k)+K2/2});

K4=h*subs(f,{t,u},{T(k)+h,Y(k)+K3));

Y(k+1)=Y(k)+(K1+2*K2+2*K3+K4)/6；

end

R[T'Y']；

h=0.1时，得到的解:

t近似解精确解

00.000.00

0.000.2800.560

0.000.470.69

0.000.2230.261

0.000.870.61

0.000.090.87

0.000.270.54

0.000.0970.329

0.000.5480.519

0.000.9920.159

1.000.8250.027

1.000.8460.480

1.000.1770.124

1.000.410.07

1.000.090.83

1.000.490.58

1.000.720.22

1.000.820.59

1.000.0500.056

1.000.010.33

2.000.640.34

h=0.05时得到的解:

t近似解精确解

00.000.00

0.000.8840.539

0.000.1210.560

0.0000.3970.043

0.000.430.69

0.000.870.42

0.000.2190.261

0.000.530.13

0.000.820.61

0.000.570.33

0.000.860.87

0.000.2350.275

0.000.210.54

0.000.600.51

0.000.2780.329

0.000.0230.617

0.000.6380.519

0.000.5120.010

0.000.7010.159

0.000.4030.377

1.000.9070.027

1.000.7040.093

1.000.3140.480

1.0000.100.50

1.000.1850.124

1.000.100.26

1.000.120.07

1.000.010.64

1.000.370.83

1.000.030.43

1.000.300.58

1.000.690.37

1.000.530.22

1.000.660.17

1.000.700.59

1.000.770.29

1.000.0590.056

1.000.780.16

1.000.660.33

1.000.0400.040

2.000.420.34

当h减半时，（1）中的全局误差缩小，和最终的预期相符。

近似解与精确解图形:

5.设

2

-1

Tn=

-1

2

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

Tn+21n—1n

—InTn+27n

~In

—InTn+2In

其中心为加介单位矩阵。设工=(1.1,1，…，1)7为线性方程组=/＞的真解,右边的

向量/曲这个真解给出。

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

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

E=IN

其中”取值为4,5,6.

解答：

(l)Cholesky分解程序代码(Ma11ab)

只列出n等于4时的代码,n等于5和6时代码类似。

c1ear

clc

%定义T(n)

T4=[2-100;-12-10;0-12-1；00-12]；

%定义单位矩阵

14=eye(4)；

/定义0矩阵

Z4=zeros(4,4);

/定义A(n)

%定义A(n)中的对角元素

d4=T4+2*14;

A4=[d4-14Z4Z4;-I4d4-14Z4;Z4-I4d4-I4；Z4Z4-

14d4];

%定义x

x4=ones(16,1)；

%求右边向量b

b4=A4*x4;

%(1)用Cho1esky分解法求解方程组

14=zeros(16,16);

%n为4时求解L

n=16;

forj=l:n

sum=0;

fork=l:j-1

sum=sum+l4(j,k)14(j,k);

end

14(j,j)=sqrt(A4(j,j)-sum);

fori=j+l:n

sum=0;

fork=l:j-1

sum=sum+14(i,k)*14(j,k);

end

14(i,j)=(A4(i,j)—sum)/14(j,j)

end

end

%n等于4时求解方程的解LL'x=b

%Ly=bL'x=y

%计算y

y4=zeros(16,1);

y4(l)=b4(1)/14(1,1);

fori=1:16

sum=0;

fork=1:i—1

sum=sum+l4(i,k)*y4(k);

end

y4(i)=(b4(i)-sum)/l4(i,i);

end

%计算x_4

x_4=zeros(16,1);

x_4(l6)=y4(l6)/14(16,16);

fori=15：-l:l

sum=0;

fork=i+l：16

sum=sum+14(k,i)*x_4(k);

end

x_4(i)=(y4(i)-sum)/14(i,i);

end

运营结果：(n等于4,5,6时)x的转置

n=4时：［1111111111111111］

n=5时：［1111111111111111111111

1111

n=6时：［111111111111111111111

111111111111111］

(2)Jacobi迭代法代码(只列出n等于4时的代码，n等于5和6时类似)

%取初始向量为0

jx40=zeros(16,1);

jx41=zeros(16,1);

dx4=zeros(16,1);

while1

fori=l:16

sum=0;

forj=l:16

sum=sum+A4(i,j)*jx40(j);

end

dx4(i)=(b4(i)—sum)/A4(i,i)；

jx41(i)=jx40(i)+dx4(i);

end

ifnorm(jx41-jx40)>le—8

%迭代

jx40=jx41;

else

break；

end

end

disp(jx41');

运营结果：(n等于4,5,6时)x的转置

当n=4时：

[0.770.460.460.770.460.230.230.460.460.230.

230.460.770.460.460.77]

当n=5时：

E0.790.500.580.500.790.500.370.000.370.5

00.580.000.160.000.580.500.370.000.370.500.79

0.500.580.500.79]

当n=6时的结果：

[0.840.430.830.830.430.840.43

0.670.260.260.670.43

0.830.260.100.100.260.830.830.26

0.100.100.260.83

0.430.670.260.260.670.430.840.4

30.830.830.430.84]

(3)Gauss-Seide1迭代法代码(只列出n等于4时的代码，n等于5和6

时类似)

%取初始向量为0

gsx40=zeros(16,1);

gsx41=zeros(16,1);

gs—dx4=zeros(16,1);

while1

fori=1:16

sum=0;

forj=l：i-1

sum=sum+A4(i,j)*gsx40(j);

end

forj=i:16

sum=sum+A4(i,j)*gsx40(j);

end

gs_dx4(i)=(b4(i)-sum)/A4(i,i)；

gsx41(i)=gsx40(i)+gs_dx4(i)；

end

ifnorm(gsx41-gsx40)>le—8

%迭代

gsx40=gsx41;

else

break;

end

end

disp(gsx41');

运营结果：（n等于4,5,6时）x的转置

当n=4时的结果:

[0.770.460.460.77

0.460.230.230.46

0.460.230.230.46

0.770.460.460.77]

当n=5时的结果:

[0.790.500.