同济大学数值分析matlab编程.docx
- 文档编号:8881804
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:18
- 大小:72.39KB
同济大学数值分析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;x=x1;
while(abs(feval(f,x))>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;
x=x0;
while(abs(feval(f,x))>tol)
x0=x;
x=x0-feval(f,x0)./feval(g,x0);
iter=iter+1;
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
其他
>>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
1牛顿法解方程组根
Functions=NewtonIterate(x,eps)
%Newton迭代法求解非线性方程组的解
%x为迭代初值,eps为允许误差
Ifnargin==1
eps=1.0e-6;
elseifnargin<1
return
end
x1=fx1(x);%非线性方程组函数fx1.m
x2=-dfx1(x);%非线性方程组的导数函数dfx1.m
x0=x2\x1';
whilenorm(x0)>=eps
x=x0'+x;
x1=fx1(x);
x2=-dfx1(x);
x0=x2\x1';
end
s=x0'+x;
return
functiony=fx1(x)%函数1
y=x^3-2*x-5
functiony=dfx1(x)%函数2
y=3*x^2-2
%命令行
>>x=NewtonIterate(2.5,eps)
x=1307/624
2拉格朗日插值法
functionyh=lagrange(x,y,xh)
n=length(x);
m=length(xh);
yh=zeros(1,m);
c1=ones(n-1,1);
c2=ones(1,m);
fori=1:
n
xp=x([1:
i-1i+1:
n]);
yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));
end
3牛顿插值法
function[p,q]=chashang(x,y)
n=length(x);
p(:
1)=x;
p(:
2)=y;
forj=3:
n+1
p(1:
n+2-j,j)=diff(p(1:
n+3-j,j-1))./(x(j-1:
n)-x(1:
n+2-j));
end
q=p(1,2:
n+1)';
4改进乘幂法
function[t,y]=eigIPower(a,xinit,ep)
v0=xinit;
[t,ti]=max(abs(v0));
lam0=v0(ti);
u0=v0/lam0;
flag=0;
while(flag==0)
v1=a*u0;
[tv,ti]=max(abs(v1));
lam1=v1(ti);
u0=v1/lam1;
err=abs(lam0-lam1);
if(err<=eo)
flag=1;
end
lam0=lam1;
end
t=lam1;
y=u0;
5高斯赛德尔迭代法求解方程组
function[x,iter]=gs(A,b,tol)
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
x=zeros(size(b));
foriter=1:
500
x=(D-L)\(b+U*x);
error=norm(b-A*x)/norm(b);
if(error break; end end 6复合梯形公式 7复合辛普森公式 8欧拉公式求解微分方程 function[x,y]=odeEuler(f,y0,a,b,n) y (1)=y0; h=(b-a)/n; x=a: h: b; fori=1: n, y(i+1)=y(i)+h*feval(f,x(i),y(i)); end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 同济大学 数值 分析 matlab 编程