工程数学线性方程组直接解法数值试验.docx
《工程数学线性方程组直接解法数值试验.docx》由会员分享，可在线阅读，更多相关《工程数学线性方程组直接解法数值试验.docx（11页珍藏版）》请在冰点文库上搜索。
工程数学线性方程组直接解法数值试验
functionx=gauss1（a,b）
n=length（b）;
fork=1:
n-1
fori=k+1:
n
m（i,k）=a（i,k）/a（k,k）;
a（i,k）=0;
b（i）=b（i）-m（i,k）*b（k）;
forj=k+1:
n
a（i,j）=a（i,j）-a（k,j）*m（i,k）;
end
end
end
x（n）=b（n）/a（n,n）;
fori=n-1:
-1:
1
s=0;
forj=i+1:
n
s=s+a（i,j）*x（j）;
end
x（i）=（b（i）-s）/a（i,i）;
end
x=x';
a=[234;352;4330];b=[6,5,32]'
a（1,2）
a=[123;456;788];b=[61523];
x=gauss1（a,b）
x=
1
1
1
x=
NaN
NaN
NaN
x=
-13
8
2
a=[2-11;-1-23;131]
b=[456]’
x=gauss1（a,b）
x=gauss1（a,b）
x=
1.11111111111111
0.777777777777778
2.55555555555556
a=[12-33;-1831;111]
b=[15-156]’
x=gauss1（A,b））
m=
0
1.5
2
a=
234
00.5-4
0-322
b=
6
-4
20
m=
00
1.50
2-6
a=
234
00.5-4
00-2
b=
6
-4
-4
x=
-13
8
2
functionx=gauss4（a,b）
a=[a,b];
n=length（b）;
fork=1:
n-1
s=a（k,k）;
p=k;
fori=k+1:
n
ifabs（s）s=a（i,k）;
p=i;
end
end
ifp~=k
forj=k:
n+1
t=a（k,j）;
a（k,j）=a（p,j）;
a（p,j）=t;
end
end
a,p
fori=k+1:
n
m（i,k）=a（i,k）/a（k,k）;
a（i,k）=0;
forj=k+1:
n+1
a（i,j）=a（i,j）-a（k,j）*m（i,k）;
end
end
end
a
x（n）=a（n,n+1）/a（n,n）;
fori=n-1:
-1:
1
s=0;
forj=i+1:
n
s=s+a（i,j）*x（j）;
end
x（i）=（a（i,n+1）-s）/a（i,i）;
end
x=x';
a=[234;352;4330];b=[6,5,32]'
x=gauss4（a,b）
a=[234;352;4330];b=[6,5,32]'
x=gauss4（a,b）
b=
6
5
32
a=
433032
3525
2346
p=
3
a=
433032
02.75-20.5-19
01.5-11-10
p=
2
a=
433032
02.75-20.5-19
000.181********81820.363636363636363
x=
-13
8
2
x=
-13
8
2
functionL=sanjiao3（a）
[n,n]=size（a）;
a（1,1）=sqrt（a（1,1））;
fori=2:
n
a（i,1）=a（i,1）/a（1,1）;
end
fork=2:
n-1
a（k,k）=a（k,k）-dot（a（k,1:
k-1）,a（k,1:
k-1））;
a（k,k）=sqrt（a（k,k））;
fori=k+1:
n
a（i,k）=（a（i,k）-dot（a（i,1:
k-1）,a（k,1:
k-1）））/a（k,k）;
end
end
a（n,n）=a（n,n）-dot（a（n,1:
n-1）,a（n,1:
n-1））;
a（n,n）=sqrt（a（n,n））;
L=tril（a）;
a=[9189-27;18450-45;901269;-27-459135]
a=
9189-27
18450-45
901269
-27-459135
>>L=sanjiao3（a）
L=
3000
6300
3-690
-9363
>>L*L'
ans=
9189-27
18450-45
901269
-27-459135
functionx=sanduijiao（a,b,c,f）
n=length（a）;
q
（1）=a
（1）;
b=[0,b];
fori=2:
n
p（i）=b（i）/q（i-1）;
q（i）=a（i）-p（i）*c（i-1）;
end
fori=2:
n
f（i）=f（i）-p（i）*f（i-1）;
end
x（n）=f（n）/q（n）;
fori=n-1:
-1:
1
x（i）=（f（i）-c（i）*x（i+1））/q（i）;
end
a=[22222];
>>b=[1111];
>>c=[-1-1-1-1];
f=[12223];x=sanduijiao（a,b,c,f）
x=
11111
function[L,U]=sanjiao1（a）
[n,n]=size（a）;
fori=2:
n
a（i,1）=a（i,1）/a（1,1）;
end
fork=2:
n
forj=k:
n
a（k,j）=a（k,j）-dot（a（k,1:
k-1）,a（1:
k-1,j））;
end
fori=k+1:
n
a（i,k）=（a（i,k）-dot（a（i,1:
k-1）,a（1:
k-1,k）））/a（k,k）;
end
end
U=triu（a）;
L=tril（a,-1）+eye（n）;
a=[9189-27;18450-45;901269;-27-459135]
function[L,U,P]=sanjiao2（a）
[n,n]=size（a）;
I=[1:
n];
fork=1:
n
forj=k:
n
a（j,k）=a（j,k）-sum1（a（j,:
）,a（:
k）,1,k-1）;
s（j）=a（j,k）;
end
m=s（k）;
t=k;
forj=k+1:
n
ifabs（m）m=s（j）;
t=j;
end
end
p=I（k）;I（k）=I（t）;I（t）=p;
forj=1:
n
p=a（k,j）;a（k,j）=a（t,j）;a（t,j）=p;
end
forj=k+1:
n
a（j,k）=a（j,k）/a（k,k）;
end
forj=k+1:
n
a（k,j）=a（k,j）-sum1（a（k,:
）,a（:
j）,1,k-1）;
end
end
L=tril（a,-1）+eye（n）;
U=triu（a）;
fori=1:
n
P（i,I（i））=1;
end
function[L,U]=sanjiao1（a）
[n,n]=size（a）;
fori=2:
n
a（i,1）=a（i,1）/a（1,1）;
end
fork=2:
n
forj=k:
n
a（k,j）=a（k,j）-dot（a（k,1:
k-1）,a（1:
k-1,j））;
end
fori=k+1:
n
a（i,k）=（a（i,k）-dot（a（i,1:
k-1）,a（1:
k-1,k）））/a（k,k）;
end
end
U=triu（a）;
L=tril（a,-1）+eye（n）;
升级会员 