大学数学实验5参考答案

实验5 线性方程组的解法

I. 实验内容及要点

一. 注意以下函数的用法

1. break：可以导致包含该命令的while、for循环终止。也可以在if-end、switch-case、

try-catch中导致中断

2.

3.

4.

5.

二.

三. continue：跳过位于其后的循环中的其他命令，执行循环的下一次迭代。 return：结束该命令所在函数的执行，把控制交给主调函数或命令窗口。 error（’message’）：显示出错信息message，终止程序。 warning(‘message’)：显示警告信息message，程序继续运行。 注意直接法中误差的判断（条件数的应用）和迭代法中收敛性的判断（见函数isshoulian） LU消元法的程序

function x=xiaoyuan(a,b)

[m n]=size(a); %可以讨论m n 的大小关系

[l u]=lu(a);

s=inv(l)*[a,b];

x=ones(m,1);

for i=m:-1:1

h=s(i,m+1);

for j=m:-1:1 %

if j~=i % h=s(i,m+1)-

h=h-x(j)*s(i,j); % (s(i,1:m) *x

end % -s(i,i))

end

x(i)=h/s(i,i);

end

四. function y=isshoulian(a)

s=size(a);if s(1)~=s(2),error('请输入方阵'),end

n=s(1);

for i=1:n

m=0;

for j=1:n

if j~=i

m=m+abs(a(i,j));

end

end

if abs(a(i,i))<m

y=0; %迭代不收敛

return

end

end

y=1;%迭代收敛

五. 雅克比迭代

function x=ydiedai(a,b,n)

if isshoulian(a)==0

warning('迭代不收敛')

return

end

l=length(b);t=b;b=zeros(l,1); %确保参与运算的是列向量

for i=1:l

b(i)=t(i);

end

[m m]=size(a);

d=diag(diag(a));

l=-tril(a,-1); %或 l=-tril(a)+d;

u=-triu(a,1); %或 u=-triu(a)+d;

b1=inv(d)*(l+u);

f1=inv(d)*b;

x=zeros(m,1);

for i=1:n %常用while循环来设计带误差的终止条件

x=b1*x+f1;

end

六. 高斯——赛德尔迭代

function x=gdiedai(A,b,x0,tol)

l1=length(x0);h=zeros(l1,1); % x0=x0(:);

for i=1:l1

h(i)=x0(i);

end

l2=length(b);t=b;b=zeros(l2,1); %b=b(:);

for i=1:l2

b(i)=t(i);

end

[m n]=size(A); %.....

d=diag(diag(A));

l=-tril(A,-1); %或 l=-tril(A)+d;

u=-triu(A,1); %或 u=-triu(A)+d;

b2=inv(d-l)*u;

f2=inv(d-l)*b;

x1=h; %即x0

x=b2*x1+f2;

i=1;

while abs(x-x1)>tol %常用范数来做判断

x1=x;i=i+1;

x=b2*x1+f2;

end

x;i

II. 课后作业

一.

二. 略 解：

1. a=[3.021 2.714 6.913;1.031 -4.273 1.121;5.084 -5.832 9.155];

b=[12.648;-2.121;8.407];

h=det(a) %判断a是否几近奇异，进而判断是否可能病态

x=xiaoyuan(a,b)

a(2,2)=-4.275;

h1=det(a) %判断a是否几近奇异，进而判断是否可能病态

x=xiaoyuan(a,b)

2.

3. 解： a=hilb(10);

x=ones(10,1);

b=a*x;

b=b.*(1+0.01);

x1=xiaoyuan(a,b);

x2=gdiedai(a,b,x,0.001);

x3=ydiedai(a,b,3);

[b x1 x2 x3]

c=cond(a)

p=max(abs(eig(a)))

三. 解：

n=1000;b=[1:n]';

a1=sparse(1:n,1:n,4,n,n);

a2=sparse(2:n,1:n-1,1,n,n);

a=a1+a2+a2';

% 输出用稀疏矩阵求解的时间t1

tic; x=a\b; t1=toc

% 与满阵做比较

aa=full(a);

% 输出用满阵求解的时间

tic; xx=aa\b; t2=toc

% 为检验x与xx是否相同分别输出其分量之和

y=sum(x)

yy=sum(xx)

四. 解：

1. % 本题可以转化为求解方程组

% 例如a题，可转化为

t1*sin(20*pi/180)=5;w+t2*sin(10*pi/180)=5;t1*cos(20*pi/180)-t1*cos(10*pi/180)=0 % 以下求解a

a=[sin(20*pi/180) 0 0;0 sin(10*pi/180) 1;cos(20*pi/180) -cos(10*pi/180) 0]; b=[5;5;0]

x1=xiaoyuan(a,b)

x2=gdiedai(a,b,[0 0 0],0.001)%看结果如何，若不行，就看范数

x3=ydiedai(a,b,3) %是否小于1，即迭代是否收敛

2. % 本题可以转化为求解方程组

% 例如b题，可转化为

l1*sin(20*pi/180)+l2*sin(10*pi/180)=d;l1*cos(20*pi/180)+l2*cos(10*pi/180)=h % 以下求解b题

a=[sin(20*pi/180),sin(10*pi/180);cos(20*pi/180),cos(10*pi/180)];

d=2;h=8;b=[d;h];

L=xiaoyuan(a,b)

3. % 本题可以转化为求解方程组

% 例如c题，可转化为

% t1*sin(40*pi/180)=5;

% t2*sin(30*pi/180)+w1=5;

% t1*cos(40*pi/180)-t2*cos(30*pi/180)=0;

% t2*sin(30*pi/180)-t3*sin(20*pi/180)-w2=0;

% t2*cos(30*pi/180)-t3*cos(20*pi/180)=0;

% 以下求解c题

a=[sin(40*pi/180) 0 0 0 0;0 sin(30*pi/180) 0 1 0;cos(40*pi/180) -cos(30*pi/180) 0 0 0; 0 sin(30*pi/180) -sin(20*pi/180) 0 -1;0 cos(30*pi/180) -cos(20*pi/180) 0 0];

b=[5;5;0;0;0];

xiaoyuan(a,b)

五.

六.

七.

八.

九.

十. 略 略 略 略 略 略

本文来源：https://www.bwwdw.com/article/kqnj.html