同济大学数值分析报告matlab编程题总汇编Word文档格式.docx
《同济大学数值分析报告matlab编程题总汇编Word文档格式.docx》由会员分享，可在线阅读，更多相关《同济大学数值分析报告matlab编程题总汇编Word文档格式.docx（19页珍藏版）》请在冰豆网上搜索。
x=x1-feval（f,x1）.*（x1-x0）./（feval（f,x1）-feval（f,x0））;
x0=x1;
x1=x;
end
editf.m
functionv=f（x）
v=x.*log（x）-1;
editg.m
functionz=g（y）
z=y.^5+y-1;
[x1iter1]=gexianfa（'
f'
1,3,1e-10）
x1=
1.7632
iter1=
6
[x2iter2]=gexianfa（'
g'
0,1,1e-10）
x2=
0.7549
iter2=
8
3.使用GS迭代求解下述线性代数方程组：
editgsdiedai.m
function[xiter]=gsdiedai（A,x0,b,tol）
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
x=x0;
while（（norm（b-A*x）./norm（b））>
tol）
x0=x;
x=（D-L）\（U*x0+b）;
A=[521;
-142;
1-310];
b=[-12103]'
tol=1e-4;
x0=[000]'
[xiter]=gsdiedai（A,x0,b,tol）;
x
x=
-3.0910
1.2372
0.9802
iter
iter=
4.用四阶Range-kutta方法求解下述常微分方程初值问题（取步长h=0.01）
editksf2.m
functionv=ksf2（x,y）
v=y+exp（x）+x.*y;
a=1;
b=2;
h=0.01;
n=（b-a）./h;
x=[1:
0.01:
2];
y
（1）=2;
fori=2:
（n+1）
k1=h*ksf2（x（i-1）,y（i-1））;
k2=h*ksf2（x（i-1）+0.5*h,y（i-1）+0.5*k1）;
k3=h*ksf2（x（i-1）+0.5*h,y（i-1）+0.5*k2）;
k4=h*ksf2（x（i-1）+h,y（i-1）+k3）;
y（i）=y（i-1）+（k1+2*k2+2*k3+k4）./6;
y
调用函数方法
editRangekutta.m
function[xy]=Rangekutta（f,a,b,h,y0）
x=[a:
h:
b];
n=（b-a）/h;
y
（1）=y0;
k1=h*（feval（f,x（i-1）,y（i-1）））;
k2=h*（feval（f,x（i-1）+0.5*h,y（i-1）+0.5*k1））;
k3=h*（feval（f,x（i-1）+0.5*h,y（i-1）+0.5*k2））;
k4=h*（feval（f,x（i-1）+h,y（i-1）+k3））;
[xy]=Rangekutta（'
ksf2'
1,2,0.01,2）;
5.取
，请编写Matlab程序，分别用欧拉方法、改进欧拉方法在
上求解初值问题。
editEuler.m
function[xy]=Euler（f,a,b,h,y0）
n=（b-a）./h;
y（i）=y（i-1）+h*feval（f,x（i-1）,y（i-1））;
editgaijinEuler.m
function[xy]=gaijinEuler（f,a,b,h,y0）
y1=y（i-1）+h*feval（f,x（i-1）,y（i-1））;
y2=y（i-1）+h*feval（f,x（i）,y1）;
y（i）=（y1+y2）./2;
editksf3.m
functionv=ksf3（x,y）
v=x.^3-y./x;
[xy]=Euler（'
ksf3'
1,2,0.2,0.4）
1.00001.20001.40001.60001.80002.0000
y=
0.40000.52000.77891.21651.88362.8407
[xy]=gaijinEuler（'
0.40000.58950.92781.46152.24643.3466
6.请编写复合梯形积分公式的Matlab程序，计算下面积分的近似值，区间等分
。
编写辛普森积分公式的Matlab程序，计算下面积分的近似值，区间等分
、
edittixingjifen.m
functions=tixingjifen（f,a,b,n）
x=linspace（a,b,（n+1））;
y=zeros（1,length（x））;
y=feval（f,x）
h=（b-a）./n;
s=0.5*h*（y
（1）+2*sum（y（2:
n））+y（n+1））;
editsimpson.m
functionI=simpson（f,a,b,n）
h=（b-a）/n;
x=linspace（a,b,2*n+1）;
y=feval（f,x）;
I=（h/6）*（y
（1）+2*sum（y（3:
2:
2*n-1））+4*sum（y（2:
2*n））+y（2*n+1））;
editksf4.m
functionv=ksf4（x）
v=1./（x.^2+1）;
tixingjifen（'
ksf4'
0,1,20）
ans=
0.7853
simpson（'
0,1,10）
0.7854
editksf5.m
functionv=ksf5（x）
if（x==0）
v=1;
else
v=sin（x）./x;
（第二个函数‘ksf5’调用求积函数时，总显示有错误：
“NaN”，还没调试好。
见谅！
7.用
迭代方法对下面方程组求解，取初始向量
editJacobi.m
function[xiter]=Jacobi（A,x0,b,tol）
while（norm（A*x-b）/norm（b）>
x=D\（（L+U）*x0+b）;
A=[24-4;
333;
442];
b=[2-3-2]'
x0=[32-1]'
[x,iter]=Jacobi（A,x0,b,1e-4）
1
-1
3
8.用牛顿法求解方程
在
附近的根。
editNewton.m
function[xiter]=Newton（f,g,x0,tol）
done=0
while~done
x=x0-feval（f,x0）/feval（g,x0）;
done=norm（x-x0）<
=tol;
if~done,x0=x;
end
editksf6.m
functionv=ksf6（x）
v=x*cos（x）+2;
editksg6.m
functionz=ksg（y）
[xiter]=Newton（'
ksf6'
'
ksg6'
2,1e-4）
2.4988
9.分别用改进乘幂法、反幂法计算矩阵A的按模最大特征值及其对应的特征向量、按模最小特征值及其对应的特征向量。
editep.m
function[t,x]=ep（A,x0,tol）
[tv0ti0]=max（abs（x0））;
lam0=x0（ti0）;
x0=x0./lam0;
x1=A*x0;
[tv1ti1]=max（abs（x1））;
lam1=x1（ti1）;
x1=x1./lam1;
while（abs（lam0-lam1）>
lam0=lam1;
[tv1ti1]=max（abs（x1））;
lam1=x1（ti1）;
t=lam1;
x=x1;
editfanep.m
function[t,x]=fanep（A,x0,tol）
x1=A\x0;
while（abs（1/lam0-1/lam1）>
t=1/lam1;
A=[126-6;
6162;
-6216];
x0=[10.5-0.5]'
[t,x]=ep（A,x0,tol）
t=
21.5440
1.0000
0.7953
-0.7953
[t,x]=fanep（A,x0,tol）
4.4560
-0.6287
0.6287
10.将积分区间n等分，用复合梯形求积公式计算定积分
比较不同
值时的误差（画出平面上log（n）-log（Error）图）．
editksf7.m
functionv=ksf7（x）
v=sqrt（1+x.^2）;
I=quad（'
ksf7'
1,3）;
n=1:
100;
fori=1:
100
x（i）=tixingjifen（'
1,3,i）;
error（i）=abs（I-x（i））;
plot（log10（n）,log10（error（n）））
11.用
editsor.m
function[x,iter]=sor（A,x0,b,omega,tol）
x=（D-omega*L）\（omega*b+（1-omega）*D*x+omega*U*x）;
A=[2-10;
-13-1;
0-12];
b=[18-5]'
x0=[10-1]'
[x,iter]=sor（A,x0,b,1.1,1e-4）
1.9999
3.0000
-1.0000
5
12.编程实现求解满足下列条件的区间[-1,2]上的三次样条函数S（x）,并画出此样条函数的图形:
xi-1012
f（xi）-1010
f（xi）'
0-1
functionspl
x=[-1012]
y=[0-1010-1]
pp=csape（x,y,'
complete'
[breaks,coefs,npolys,ncoefs,dim]=unmkpp（pp）
2
if-1<
=xh<
=0
yh=coefs（1,1）*（xh+1）.^3+coefs（1,2）*（xh+1）.^2+coefs（1,3）*（xh+1）+coefs（1,4）
elseif0<
=1
yh=coefs（2,1）*（xh）.^3+coefs（2,2）*（xh）.^2+coefs（2,3）*（xh）+coefs（2,4）
elseif1<
=2
yh=coefs（3,1）*（xh-1）.^3+coefs（3,2）*（xh-1）.^2+coefs（3,3）*（xh-1）+coefs（3,4）
return
plot（xh,yh,'
13.二分法程序
functionx=bisect（f,a,b,tol）
ifnargin<
4
tol=1e-12
fa=feval（f,a）
fb=feval（f,b）
whileabs（a-b）>
tol
x=（a+b）/2
fx=feval（f,x）
ifsign（fx）==sign（fa）
a=x
fa=fx
elseifsign（fx）==sign（fb）
b=x
fb=fx
else
其他
A=[1234;
5678];
A
A=
1234
5678
[mn]=size（A）
m=
2
n=
4
x=[12345];
12345
length（x）
A=[234;
119;
12-6];
234
119
12-6
[LU]=lu（A）;
L
L=
1.000000
0.50001.00000
0.5000-1.00001.0000
U
U=
2.00003.00004.0000
0-0.50007.0000
00-1.0000
x=linspace（0,1,11）;
x'
0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
x=[1234];
y=[6111827];
p=polyfit（x,y,2）
p=
1.00002.00003.0000
diag（ones（4,1）,1）
01000
00100
00010
00001
00000
=1
升级会员 