线性方程组与MATLAB应用.docx
《线性方程组与MATLAB应用.docx》由会员分享，可在线阅读，更多相关《线性方程组与MATLAB应用.docx（12页珍藏版）》请在冰豆网上搜索。
线性方程组与MATLAB应用
实验5线性方程组与MATLAB应用
一、实验名称：
线性方程组与MATLAB应用。
二、实验目的：
理解线性方程组直接法与迭代法思想，掌握常用算法的设计，掌握用MATLAB实现的数值解法。
三、实验要求：
1、编写列主元消去法程序，并举例子。
编写ＬＵ分解法程序，并举例子。
对两种算法作出对比。
利用MATLAB函数“\”计算上面例子。
2、编写Jacobi或Gauss-Seidel迭代法程序，并举例子。
编写SOR迭代法程序，并举例子选择较佳的松弛因子。
对比说明两种迭代法的收敛速度。
3、谈一谈你对线性方程组直接法与迭代法的理解。
四、实验内容：
1、利用列主元法解方程：
gaussh.m文件为：
functionx=gauss（a,b）
%%判断是否可以求出解
temp1=size（a）;temp2=size（b）;
iftemp1
（1）~=temp1
（2）|temp1
（1）~=temp2
（1）|temp2
（2）~=1|det（a）==0
error（'矩阵或向量的大小不对应或矩阵为奇异矩阵'）
return;
end
%%化为三角矩阵过程
leng=temp1
（1）;
fori=1:
leng-1
pos=i;num=abs（a（i,i））;
forj=i+1:
leng
ifabs（a（j,i））>num
num=abs（a（j,i））;
pos=j;
end
end
ifpos~=i
fork=i:
leng
temp=a（i,k）;
a（i,k）=a（pos,k）;
a（pos,k）=temp;
end
temp=b（pos）;
b（pos）=b（i）;
b（i）=temp;
end
forj=i+1:
leng
temp=a（j,i）/a（i,i）;
forh=i:
leng
a（j,h）=a（j,h）-a（i,h）*temp;
end
b（j）=b（j）-temp*b（i）;
end
end
%%进行回代
x（leng）=b（leng）/a（leng,leng）;
ifleng>1
fori=leng-1:
-1:
1
x（i）=b（i）;
forj=i+1:
leng;
x（i）=x（i）-a（i,j）*x（j）;
end
x（i）=x（i）/a（i,i）;
end
end
x=x';
控制命令为：
>>clear
>>a=[123;252;315];b=[141820]';
>>gaussh（a,b）
ans=
1.0000
2.0000
3.0000
利用LU列主元法解方程
mylu.m文件内容为：
function[l,u,p]=mylu（a）
siz=size（a）;
ifsiz
（1）~=siz
（2）
error（'矩阵非为方阵'）
return
end
ifdet（a）==0
error（'矩阵不能被三角分解'）
end
si=siz
（1）;
u=a;p=eye（si）;l=eye（si）;
fori=1:
si
forj=i:
si
t（j）=u（j,i）;
fork=1:
i-1
t（j）=t（j）-u（j,k）*u（k,i）;
end
end
pos=i;ma=abs（t（i））;
forj=i+1:
si
ifmama=abs（t（j））;
pos=j;
end
end
ifpos~=i
forj=1:
si
temp=u（i,j）;
u（i,j）=u（pos,j）;
u（pos,j）=temp;
end
forj=1:
si
temp=p（i,j）;
p（i,j）=p（pos,j）;
p（pos,j）=temp;
end
temp=t（pos）;
t（pos）=t（i）;
t（i）=temp;
end
u（i,i）=t（i）;
forj=i+1:
si
u（j,i）=t（j）/t（i）;
end
forj=i+1:
si
fork=1:
i-1
u（i,j）=u（i,j）-u（i,k）*u（k,j）;
end
end
end
l=tril（u,-1）+eye（si）;
u=triu（u,0）;
控制命令为：
a=[123;252;315];b=[141820]';
>>[l,u,p]=mylu（a）
l=
1.000000
0.66671.00000
0.33330.38461.0000
u=
3.00001.00005.0000
04.3333-1.3333
001.8462
p=
001
010
100
>>p*b
ans=
20
18
14
y=l\ans
y=
20.0000
4.6667
5.5385
>>x=u\y
x=
1.0000
2.0000
3.0000
利用‘\’直接求解为：
>x=a\b
x=
1.0000
2.0000
3.0000
2、用Jacobi解方程：
…………①
Jacobi.m文件为：
function[x,n]=jacobi（a,b,E）
temp1=size（a）;temp2=size（b）;
iftemp1
（1）~=temp1
（2）|temp1
（1）~=temp2
（1）|temp2
（2）~=1|det（a）==0
error（'矩阵或向量的大小不对应或矩阵为奇异矩阵'）
return;
end
%%化为三角矩阵过程
leng=temp1
（1）;N=-tril（a,-1）-triu（a,1）;M=a+N;
t=1;
fori=1:
leng
t=t*M（i,i）;
end
ift==0
error（'对角元素都为零！
'）
end
x=zeros（leng,1）;
t=0;x=zeros（leng,1）;
x2=N*x+b;
fori=1:
leng
x1=x2（i）/M（i,i）;
t=t+（x1-x（i））^2;
x（i）=x1;
end
e=E^2;
n=1;
whilet>e
t=0;
x2=N*x+b;
fori=1:
leng
x1=x2（i）/M（i,i）;
t=t+（x1-x（i））^2;
x（i）=x1;
end
n=n+1;
ift>1000
break;
error（'不收敛'）
end
end
控制窗口命令为：
a=[8-32;411-1;6312];b=[203336]';
>>[x,n]=jacobi（a,b,0.001）
x=
3.0003
1.9999
0.9997
n=
9
用Gauss-Seidel迭代法解方程①
Gauss_seidel.m文件内容为：
function[x,n]=gauss_seidel（a,b,E）
temp1=size（a）;temp2=size（b）;
iftemp1
（1）~=temp1
（2）|temp1
（1）~=temp2
（1）|temp2
（2）~=1|det（a）==0
error（'矩阵或向量的大小不对应或矩阵为奇异矩阵'）
return;
end
%%化为三角矩阵过程
leng=temp1
（1）;N=-tril（a,-1）-triu（a,1）;M=a+N;
t=1;
fori=1:
leng
t=t*M（i,i）;
end
ift==0
error（'对角元素都为零！
'）
end
x=zeros（leng,1）;
t=0;x1=zeros（leng,1）;
fori=1:
leng
forj=1:
leng
ifjx（i）=N（i,j）*x（j）+x（i）;
end
ifj>i
x（i）=N（i,j）*x1（j）+x（i）;
end
end
x（i）=（x（i）+b（i））/M（i,i）;
end
fori=1:
leng
t=t+（x1（i）-x（i））^2;
end
e=E^2;
n=1;
whilet>e
t=0;
x1=x;
x=zeros（leng,1）;
fori=1:
leng
forj=1:
leng
ifjx（i）=N（i,j）*x（j）+x（i）;
end
ifj>i
x（i）=N（i,j）*x1（j）+x（i）;
end
end
x（i）=（x（i）+b（i））/M（i,i）;
t=t+（x1（i）-x（i））^2;
end
n=n+1;
ift>1000
break;
error（'不收敛'）
end
end
控制命令窗口命令符为：
a=[8-32;411-1;6312];b=[203336]';
>>[x,n]=gauss_seidel（a,b,0.001）
x=
2.9998
2.0001
1.0001
n=
5
用SOR迭代法解法方程①
SOR.m文件为
function[x,n]=SOR（a,b,E,p）
temp1=size（a）;temp2=size（b）;
iftemp1
（1）~=temp1
（2）|temp1
（1）~=temp2
（1）|temp2
（2）~=1|det（a）==0
error（'矩阵或向量的大小不对应或矩阵为奇异矩阵'）
return;
end
%%化为三角矩阵过程
leng=temp1
（1）;N=-tril（a,-1）-triu（a,1）;M=a+N;N=-a;
t=1;
fori=1:
leng
t=t*M（i,i）;
end
ift==0
error（'对角元素都为零！
'）
end
x=zeros（leng,1）;
t=0;x1=zeros（leng,1）;
fori=1:
leng
forj=1:
leng
ifjx（i）=N（i,j）*x（j）+x（i）;
end
ifj>=i
x（i）=N（i,j）*x1（j）+x（i）;
end
end
x（i）=x1（i）+（x（i）+b（i））/M（i,i）*p;
end
fori=1:
leng
t=t+（x1（i）-x（i））^2;
end
e=E^2;
n=1;
whilet>e
t=0;
x1=x;
x=zeros（leng,1）;
fori=1:
leng
temp=0;
forj=1:
leng
ifjx（i）=N（i,j）*x（j）+x（i）;
end
ifj>=i
x（i）=N（i,j）*x1（j）+x（i）;
end
end
x（i）=x1（i）+（x（i）+b（i））/M（i,i）*p;
t=t+（x1（i）-x（i））^2;
end
n=n+1;
ift>1000
break;
error（'不收敛'）
end
end
控制命令窗口命令为：
当松弛因子p=1.3时
[x,n]=SOR（a,b,0.001,1.3）
x=
升级会员 