同济大学数值分析报告matlab编程题总汇编.docx
《同济大学数值分析报告matlab编程题总汇编.docx》由会员分享，可在线阅读，更多相关《同济大学数值分析报告matlab编程题总汇编.docx（18页珍藏版）》请在冰豆网上搜索。
同济大学数值分析报告matlab编程题总汇编
MATLAB编程题库
1.下面的数据表近似地满足函数
，请适当变换成为线性最小二乘问题，编程求最好的系数
，并在同一个图上画出所有数据和函数图像.
解：
x=[-0.931-0.586-0.362-0.2130.0080.5440.6280.995]';
y=[0.3560.6060.6870.8020.8230.8010.7180.625]';
A=[xones（8,1）-x.^2.*y];
z=A\y;
a=z
（1）;b=z
（2）;c=z（3）;
xh=-1:
0.1:
1;
yh=（a.*xh+b）./（1+c.*xh.^2）;
plot（x,y,'r+',xh,yh,'b*'）
2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数精度为
的近似根，并写出调用方式：
文件一
文件二
functionv=f（x）
v=x.*log（x）-1;
functionz=g（y）
z=y.^5+y-1;
解：
>>editgexianfa.m
function[xiter]=gexianfa（f,x0,x1,tol）
iter=0;
while（norm（x1-x0）>tol）
iter=iter+1;
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）;
iter=0;
x=x0;
while（（norm（b-A*x）./norm（b））>tol）
iter=iter+1;
x0=x;
x=（D-L）\（U*x0+b）;
end
>>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=
6
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;
end
>>y
调用函数方法
>>editRangekutta.m
function[xy]=Rangekutta（f,a,b,h,y0）
x=[a:
h:
b];
n=（b-a）/h;
y
（1）=y0;
fori=2:
（n+1）
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））;
y（i）=y（i-1）+（k1+2*k2+2*k3+k4）./6;
end
>>[xy]=Rangekutta（'ksf2',1,2,0.01,2）;
>>y
5.取
，请编写Matlab程序，分别用欧拉方法、改进欧拉方法在
上求解初值问题。
解：
>>editEuler.m
function[xy]=Euler（f,a,b,h,y0）
x=[a:
h:
b];
n=（b-a）./h;
y
（1）=y0;
fori=2:
（n+1）
y（i）=y（i-1）+h*feval（f,x（i-1）,y（i-1））;
end
>>editgaijinEuler.m
function[xy]=gaijinEuler（f,a,b,h,y0）
x=[a:
h:
b];
n=（b-a）./h;
y
（1）=y0;
fori=2:
（n+1）
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;
end
>>editksf3.m
functionv=ksf3（x,y）
v=x.^3-y./x;
>>[xy]=Euler（'ksf3',1,2,0.2,0.4）
x=
1.00001.20001.40001.60001.80002.0000
y=
0.40000.52000.77891.21651.88362.8407
>>[xy]=gaijinEuler（'ksf3',1,2,0.2,0.4）
x=
1.00001.20001.40001.60001.80002.0000
y=
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））;
end
>>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:
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（'ksf4',0,1,10）
ans=
0.7854
>>editksf5.m
functionv=ksf5（x）
if（x==0）
v=1;
else
v=sin（x）./x;
end
（第二个函数‘ksf5’调用求积函数时，总显示有错误：
“NaN”，还没调试好。
见谅！
）
7.用
迭代方法对下面方程组求解，取初始向量
。
解：
>>editJacobi.m
function[xiter]=Jacobi（A,x0,b,tol）
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
x=x0;
iter=0;
while（norm（A*x-b）/norm（b）>tol）
iter=iter+1;
x0=x;
x=D\（（L+U）*x0+b）;
end
>>A=[24-4;333;442];
>>b=[2-3-2]';
>>x0=[32-1]';
>>[x,iter]=Jacobi（A,x0,b,1e-4）
x=
1
-1
-1
iter=
3
8.用牛顿法求解方程
在
附近的根。
解：
>>editNewton.m
function[xiter]=Newton（f,g,x0,tol）
iter=0;
done=0
while~done
x=x0-feval（f,x0）/feval（g,x0）;
done=norm（x-x0）<=tol;
iter=iter+1;
if~done,x0=x;end
end
>>editksf6.m
functionv=ksf6（x）
v=x*cos（x）+2;
>>editksg6.m
functionz=ksg（y）
z=y.^5+y-1;
>>[xiter]=Newton（'ksf6','ksg6',2,1e-4）
x=
2.4988
iter=
3
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）>tol）
x0=x1;
lam0=lam1;
x1=A*x0;
[tv1ti1]=max（abs（x1））;
lam1=x1（ti1）;
x1=x1./lam1;
end
t=lam1;
x=x1;
>>editfanep.m
function[t,x]=fanep（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（1/lam0-1/lam1）>tol）
x0=x1;
lam0=lam1;
x1=A\x0;
[tv1ti1]=max（abs（x1））;
lam1=x1（ti1）;
x1=x1./lam1;
end
t=1/lam1;
x=x1;
>>A=[126-6;6162;-6216];
>>x0=[10.5-0.5]';
>>tol=1e-4;
>>[t,x]=ep（A,x0,tol）
t=
21.5440
x=
1.0000
0.7953
-0.7953
>>A=[126-6;6162;-6216];
>>x0=[10.5-0.5]';
>>tol=1e-4;
[t,x]=fanep（A,x0,tol）
t=
4.4560
x=
1.0000
-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（'ksf7',1,3,i）;
error（i）=abs（I-x（i））;
end
>>plot（log10（n）,log10（error（n）））
11.用
迭代方法对下面方程组求解，取初始向量
。
解：
>>editsor.m
function[x,iter]=sor（A,x0,b,omega,tol）
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
x=x0;
iter=0;
while（norm（A*x-b）/norm（b）>tol）
iter=iter+1;
x0=x;
x=（D-omega*L）\（omega*b+（1-omega）*D*x+omega*U*x）;
end
>>A=[2-10;-13-1;0-12];
>>b=[18-5]';
>>x0=[10-1]';
>>[x,iter]=sor（A,x0,b,1.1,1e-4）
x=
1.9999
3.0000
-1.0000
iter=
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）
xh=-1:
0.1:
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<=xh<=1
yh=coefs（2,1）*（xh）.^3+coefs（2,2）*（xh）.^2+coefs（2,3）*（xh）+coefs（2,4）
elseif1<=xh<=2
yh=coefs（3,1）*（xh-1）.^3+coefs（3,2）*（xh-1）.^2+coefs（3,3）*（xh-1）+coefs（3,4）
else
return
end
plot（xh,yh,'r+'）
13.二分法程序
functionx=bisect（f,a,b,tol）
ifnargin<4
tol=1e-12
end
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
return
end
end
其他
>>A=[1234;5678];
>>A
A=
1234
5678
>>[mn]=size（A）
m=
2
n=
4
>>x=[12345];
>>x
x=
12345
>>length（x）
ans=
5
>>A=[234;119;12-6];
>>A
A=
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'
ans=
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）
ans=
01000
00100
00010
00001
00000
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）
xh=-1:
0.1:
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<=xh<=1
yh=coefs（2,1）*（xh）.^3+coefs（2,2）*（xh）.^2+coefs（2,3）*（xh）+coefs（2,4）
elseif1<=xh<=2
yh=coefs（3,1）*（xh-1）.^3+coefs（3,2）*（xh-1）.^2+coefs（3,3）*（xh-1）+coefs（3,4）
else
return
end
plot（xh,yh,'r+'）
13.二分法程序
functionx=bisect（f,a,b,tol）
ifnargin<4
tol=1e-12
end
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
return
end
end
升级会员 