计算机编程上机实验Word下载.docx
《计算机编程上机实验Word下载.docx》由会员分享，可在线阅读，更多相关《计算机编程上机实验Word下载.docx（24页珍藏版）》请在冰点文库上搜索。
x1（:
1）=x;
10000
x1（:
i）=a1^-1*a2*x1（:
i-1）+a1^-1*b;
ifabs（max（x1（:
i））-max（x1（:
i-1）））<
=c;
disp（x1）;
jacobi=x1（:
i）;
break
end
2、Gauss‐Seidel迭代法
functiongauss=gauss（a,b,c,x）
n=length（a）;
d（i,j）=a（i,j）;
elseifi<
j
u（i,j）=-a（i,j）;
elseifi>
j
l（i,j）=-a（i,j）;
u=[u;
zeros（1,n）];
l=[l,zeros（n,1）];
i）=（d-l）^-1*u*x1（:
i-1）+（d-l）^-1*b;
gauss=x1（:
在matlab主程序中需要输入系数矩阵a和b，迭代精度为c，x为计算初值。
通过输入
a=
-2100
1-210
01-21
001-2
b=[-1;
0;
0]
c=1.0000e-006
x1=[0;
x2=[0.5;
0.5;
0.5]
然后打入jacobi（a,b,c,x1）得到如下结果
00.50000.5000
000.2500
000
ans=
0.5000
0.2500
0
打入gauss（a,b,c,x1）得到如下结果
Columns1through11
0.50000.75000.81250.84380.83200.82130.81400.80910.80600.80390.8026
0.50000.62500.68750.66410.64260.62790.61830.61200.60780.60510.6034
0.50000.56250.48440.45310.43460.42260.41480.40970.40630.40410.4027
0.50000.28130.24220.22660.21730.21130.20740.20480.20320.20210.2014
Columns12through22
0.80170.80110.80070.80050.80030.80020.80010.80010.80010.80000.8000
0.60220.60140.60090.60060.60040.60030.60020.60010.60010.60000.6000
0.40180.40120.40080.40050.40030.40020.40010.40010.40010.40000.4000
0.20090.20060.20040.20020.20020.20010.20010.20000.20000.20000.2000
Columns23through29
0.80000.80000.80000.80000.80000.80000.8000
0.60000.60000.60000.60000.60000.60000.6000
0.40000.40000.40000.40000.40000.40000.4000
0.20000.20000.20000.20000.20000.20000.2000
0.8000
0.6000
0.4000
0.2000
然后再打入jacobi（a,b,c,x2）得到
0.50000.75000.7500
0.50000.50000.6250
0.50000.50000.3750
0.50000.25000.2500
0.7500
0.6250
0.3750
0.2500
最后打入gauss（a,b,c,x2）得到如下结果：
2.用Gauss列主元消去法、QR方法求解如下方程组：
一、通过matlab编程
1、Gauss列主元消去法程序如下：
functiongausslzy
a=[1,2,1,2;
2,5,3,-2;
-2,-2,3,5;
1,3,2,3];
b=[4;
7;
-1;
0];
an=[a,b];
forj=1:
fori=j:
ifabs（an（i,j））==max（abs（an（j:
n,j）））
c=an（i,:
）;
d=an（j,:
an（j,:
）=c;
an（i,:
）=d;
fori=1:
ifi>
）=an（i,:
）-an（j,:
）*（an（i,j）/an（j,j））;
x=an（1:
n,1:
n）^-1*an（:
n+1）;
disp（an）;
disp（x）;
2、QR方法程序如下：
functionqr
b=0;
b2=0;
b3=0;
b=b+a（i,1）^2;
ifi==1
ui（i）=1;
else
ui（i）=0;
b=b^0.5;
ui=ui'
;
w=a（:
1）-b*ui;
q1=dan（n）-2/（w'
*w）*w*w'
a1=q1*a;
a2=a1;
a2（:
1）=[];
a2（1,:
）=[];
n-1
b2=b2+a2（i,1）^2;
u2（i）=1;
u2（i）=0;
b2=b2^0.5;
u2=u2'
w=a2（:
1）-b2*u2;
q2=dan（n-1）-2/（w'
a3=q2*a2;
a4=a3;
a4（:
a4（1,:
b3=b3+a4（i,1）^2;
u3（i）=1;
u3（i）=0;
b3=b3^0.5;
u3=u3'
w=a4（:
1）-b3*u3;
q3=dan（n-2）-2/（w'
a5=q3*a4;
qq2=dan（n）;
forj=2:
qq2（i,j）=q2（i-1,j-1）;
a1（i,j）=a3（i-1,j-1）;
qq3=dan（n）;
forj=3:
qq3（i,j）=q3（i-2,j-2）;
a1（i,j）=a5（i-2,j-2）;
q=（qq3*qq2*q1）'
x=a1^-1*q^-1*[4;
Q=q
R=a1
x
functiondan=dan（n）
dan（i,j）=1;
三、非线性方程的迭代解法
1．
通过matlab编程，程序如下：
functionnewton=newton（x）
x
（1）=x;
1000
x（i）=x（i-1）-f（x（i-1））/f1（x（i-1））;
ifabs（x（i）-x（i-1））<
10^-6
newton=x（i）;
functionf=f（x）
f=exp（x）+2^（-x）+2*cos（x）-6;
functionf1=f1（x）
f1=exp（x）-log
（2）/2^x-2*sin（x）;
代入初值x1=0，x2=-10。
得到结果如下
x=newton（0）
x=
1.8294
x=newton（-10）
-2.9865
2．利用Newton迭代法求多项式
的所有实零点，注意重根的问题。
当求解1.1附近的三重根时候有：
x（i）=x（i-1）-3*f（x（i-1））/f1（x（i-1））;
f=x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951;
f1=4*x^3-5.4*3*x^2+10.56*2*x-8.954;
结果为：
newton（0）
1.1000
当求解2.1附近的解时候程序有：
newton
（2）
2.1000
四、数值积分
一、复化梯形公式程序如下：
functiontixing=tixing
b=8;
a=2;
forn=2:
100000
t（n）=（b-a）/2/n*f（a）;
at=（b-a）/n;
t（n）=t（n）+（b-a）/n*（f（a+j*at））;
t（n）=t（n）+（b-a）/2/n*f（b）;
ifabs（t（n）-t（n-1））<
=10^-5
tixing=t;
disp（n）;
f=1/x;
结果为:
迭代次数:
53
ans=1.3865
一、Romberg公式
functionromberg=romberg
t=zeros（20）;
t（1,1）=（b-a）/2*（f（a）+f（b））;
20
t（i,1）=0.5*t（i-1）;
forj=0:
（2^（i-1）-1）
t（i,1）=t（i,1）+（b-a）/2^i*f（a+（2*j+1）*（b-a）/2^i）;
fork=2:
ifk>
=i
t（k,i）=（4^（i-1）*t（k,i-1）-t（k-1,i-1））/（4^（i-1）-1）;
ifabs（t（i,i）-t（i+1,i+1））<
romberg=t（i+1,i+1）;
disp（i+1）;
迭代次数为16
ans=1.3863
五、插值与逼近
1．给定[−1,1]上的函数
，请做如下的插值逼近：
⑴构造等距节点分别取n=5，n=8，n=10的Lagrange插值多项式；
⑵构造分段线性取n=10的Lagrange插值多项式；
⑶取Chebyshev多项式的零点：
作插值节点构造n=10的插值多项式
f（x）和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。
（1）通过matlab编程可知
functionchazhi=chazhi（n）
symsx
a（i）=-1+2/（n-1）*（i-1）;
b（i）=1/（1+25*a（i）^2）;
l（i）=1;
g（i）=sym（'
1'
ifi~=j
l（i）=l（i）/（a（i）-a（j））;
g（i）=（x-a（j））*g（i）;
la（i）=l（i）*b（i）*g（i）;
chazhi=sum（la）;
其图见下面
（2）、matlab程序如下
functionchazhi1=chazhi1
10
a（i）=-1+2/（10-1）*（i-1）;
9
la（i）=b（i）*（x-a（i+1））/（a（i）-a（i+1））+b（i+1）*（x-a（i））/（a（i+1）-a（i））;
chazhi1=la;
xx（i-1,:
）=a（i-1）:
2/81:
a（i）;
y（i,j）=b（i）*（xx（i,j）-a（i+1））/（a（i）-a（i+1））+b（i+1）*（xx（i,j）-a（i））/（a（i+1）-a（i））;
a=-1:
0.01:
1;
yy=lag（a）;
plot（xx（1,:
）,y（1,:
）,xx（2,:
）,y（2,:
）,xx（3,:
）,y（3,:
）,xx（4,:
）,y（4,:
）,xx（5,:
）,y（5,:
）,xx（6,:
）,y（6,:
）,xx（7,:
）,y（7,:
）,xx（8,:
）,y（8,:
）,xx（9,:
）,y（9,:
）,a,yy）;
结果如图所示：
（3）、matlab程序为：
functionchazhi3=chazhi3
a（i）=cos（（（2*i-1）/20）*pi）;
chazhi3=sum（la）;
2．已知函数值
xk
1
2
3
4
5
6
7
8
yk
2.51
3.30
4.04
4.70
5.22
5.54
5.78
5.40
5.57
5.70
5.80
和边界条件：
s′（0）=0.8,s′（10）=0.2．求三次样条插值函数y=s（x）并画出其图形．
functionyangtiaochazhi=yangtiaochazhi
xk=0:
1:
10;
yk=[2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80];
h=xk
（2）-xk
（1）;
l（i）=0.5;
u（i）=0.5;
g（i）=3*（u（i）*（yk（i+1）-yk（i））/h+l（i）*（yk（i）-yk（i-1））/h）;
m0=0.8;
m10=0.2;
aa（i,j）=2;
elseifi==j+1
aa（i,j）=0.5;
elseifi==j-1
g
（1）=[];
g
（1）=g
（1）-0.5*m0;
g（9）=g（9）-0.5*m10;
m=aa^-1*g'
m=[m0,m'
];
m=[m,m10];
s=[];
s=[s;
（h+2*（x-xk（i）））/h^3*（x-xk（i+1））^2*yk（i）+（h-2*（x-xk（i+1）））/h^3*（x-xk（i））^2*yk（i+1）+（x-xk（i））*（x-xk（i+1））^2/h^2*m（i）+（x-xk（i））^2*（x-xk（i+1））/h^2*m（i+1）];
yangtiaochazhi=s;
（4*x*（x-1）^2）/5+（6948927792327497*x^2*（x-1））/9007199254740992-（33*x^2*（2*x-3））/10+（251*（2*x+1）*（x-1）^2）/100
（6341574006158381*（x-1）^2*（x-2））/9007199254740992+（6948927792327497*（x-1）*（x-2）^2）/9007199254740992+（33*（2*x-1）*（x-2）^2）/10-（101*（2*x-5）*（x-1）^2）/25
（2757506526475575*（x-2）^2*（x-3））/4503599627370496+（6341574006158381*（x-2）*（x-3）^2）/9007199254740992+（101*（2*x-3）*（x-3）^2）/25-（47*（2*x-7）*（x-2）^2）/10
（1741929571910061*（x-3）^2*（x-4））/4503599627370496+（2757506526475575*（x-3）*（x-4）^2）/4503599627370496+（47*（2*x-5）*（x-4）^2）/10-（261*（2*x-9）*（x-3）^2）/50
（6495384987431321*（x-4）^2*（x-5））/18014398509481984+（1741929571910061*（x-4）*（x-5）^2）/4503599627370496+（261*（2*x-7）*（x-5）^2）/50-（277*（2*x
升级会员 