东南大学数值分析上机题作业MATLAB版.docx
- 文档编号:4509828
- 上传时间:2022-12-01
- 格式:DOCX
- 页数:27
- 大小:335.35KB
东南大学数值分析上机题作业MATLAB版.docx
《东南大学数值分析上机题作业MATLAB版.docx》由会员分享,可在线阅读,更多相关《东南大学数值分析上机题作业MATLAB版.docx(27页珍藏版)》请在冰豆网上搜索。
东南大学数值分析上机题作业MATLAB版
11
1
题目
设
,其精确值为
。
(1)编制按从大到小的顺序
,计算SN的通用程序。
(2)编制按从小到大的顺序
,计算SN的通用程序。
(3)按两种顺序分别计算
并指出有效位数。
(编制程序时用单精度)
(4)通过本次上机题,你明白了什么
clear;
N=input('请输入N值:
');
Ac=single((3/2-1/N-1/(N+1))/2);
Snl2s=single(0);
Sns2l=single(0);
fori=2:
N
Snl2s=Snl2s+1/(i*i-1);
end
fori=N:
-1:
2
Sns2l=Sns2l+1/(i*i-1);
end
fprintf('精确值为:
%f\n',Ac);
fprintf('从大到小的顺序累加得SN=%f\n',Snl2s);
fprintf('从小到大的顺序累加得SN=%f\n',Sns2l);
disp('========================================================');
程序
>>P20T17
请输入N值:
10^2
精确值为:
从大到小的顺序累加得SN=
从小到大的顺序累加得SN=
============================================================
>>P20T17
请输入N值:
10^4
精确值为:
从大到小的顺序累加得SN=
运行结果
从小到大的顺序累加得SN=
============================================================
>>P20T17
请输入N值:
10^6
精确值为:
从大到小的顺序累加得SN=
从小到大的顺序累加得SN=
============================================================
结果分析
按从大到小的顺序,有效位数分别为:
6,4,3。
按从小到大的顺序,有效位数分别为:
5,6,6。
可以看出,不同的算法造成的误差限是不同的,好的算法可以让结果更加精确。
当采用从大到小的顺序累加的算法时,误差限随着N的增大而增大,可见在累加的过程中,误差在放大,造成结果的误差较大。
因此,采取从小到大的顺序累加得到的结果更加精确。
2
题目
(1)给定初值
及容许误差
,编制牛顿法解方程f(x)=0的通用程序。
(2)给定方程
易知其有三个根
由牛顿方法的局部收敛性可知存在
当
时,Newton迭代序列收敛于根x2*。
试确定尽可能大的
。
试取若干初始值,观察当
时Newton序列的收敛性以及收敛于哪一个根。
(3)通过本上机题,你明白了什么
程序
函数m文件:
functionFu=fu(x)
Fu=x^3/3-x;
end
函数m文件:
functionFu=dfu(x)
Fu=x^2-1;
end
用Newton法求根的通用程序
clear;
x0=input('请输入初值x0:
');
ep=input('请输入容许误差:
');
flag=1;
whileflag==1
x1=x0-fu(x0)/dfu(x0);
ifabs(x1-x0) flag=0; end x0=x1; end fprintf('方程的一个近似解为: %f\n',x0); 寻找最大δ值的程序: clear eps=input('请输入搜索精度: '); ep=input('请输入容许误差: '); flag=1; k=0; x0=0; whileflag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; whileflag1==1&&m<=10^3 x1=x0-fu(x0)/dfu(x0); ifabs(x1-x0) flag1=0; end m=m+1; x0=x1; end ifflag1==1||abs(x0)>=ep flag=0; end end fprintf('最大的sigma值为: %f\n',sigma); 运行结果 (1)寻找最大的 值。 算法为: 将初值x0在从0开始不断累加搜索精度eps,带入Newton迭代公式,直到求得的根不再收敛于0为止,此时的x0值即为最大的sigma值。 运行,得到在不同的搜索精度下的最大sigma值。 >>Find 请输入搜索精度: 10^-6 请输入容许误差: 10^-6 最大的sigma值为: >>Find 请输入搜索精度: 10^-4 请输入容许误差: 10^-6 最大的sigma值为: >>Find 请输入搜索精度: 10^-2 请输入容许误差: 10^-6 最大的sigma值为: (2)运行 在 内取初值,运行结果如下: X0 Xk -1000 -500 -100 -10 -5 可见,在 区间内取初值,Newton序列收敛,且收敛于根 。 在 内取初值,运行结果如下: X0 Xk 可见,在 内取初值,Newton序列收敛,且收敛于根 。 在 内内取初值,运行结果如下: X0 Xk 可见,在 内取初值,Newton序列收敛,且收敛于根0。 在 内取初值,运行结果如下: X0 Xk 可见,在 内取初值,Newton序列收敛,且收敛于根 在 内取初值,运行结果如下: X0 Xk 5 10 100 500 1000 可见,在 内取初值,Newton序列收敛,且收敛于根 3 题目 对于某电路的分析,归结为求解线性方程组RI=V,其中 (1)编制解n阶线性方程组 的列主元高斯消去法的通用程序; (2)用所编程序线性方程组 ,并打印出解向量,保留5位有效数字; (3)本题编程之中,你提高了哪些编程能力 程序 n=input('请输入线性方程组阶数: n='); b=zeros(1,n); A=input('请输入系数矩阵: A=\n'); b(1,: )=input('请输入线性方程组右端向量: b=\n'); b=b'; C=[A,b]; fori=1: n-1 [maximum,index]=max(abs(C(i: n,i))); index=index+i-1; T=C(index,: ); C(index,: )=C(i,: ); C(i,: )=T; fork=i+1: n ifC(k,i)~=0 C(k,: )=C(k,: )-C(k,i)/C(i,i)*C(i,: ); end end end %%回代求解 x=zeros(n,1); x(n)=C(n,n+1)/C(n,n); fori=n-1: -1: 1 x(i)=(C(i,n+1)-C(i,i+1: n)*x(i+1: n,1))/C(i,i); end disp('方程组的解为: '); fprintf('%.5g\n',x); 运行结果 运行程序,输入系数矩阵和方程组右端列向量。 运行过程与结果如下图所示: >>P126T39 请输入线性方程组阶数: n=4 请输入系数矩阵: A= [00;0;0;00] 请输入线性方程组右端向量: b= [] 方程组的解为: 2495 >>P126T39 请输入线性方程组阶数: n=9 请输入系数矩阵: A= [31-13000-10000;-1335-90-110000;0-931-1000000;00-1079-30000-9;000-3057-70-50;0000-747-3000;00000-304100;0000-50027-2;000-9000-229] 请输入线性方程组右端向量: b= [-1527-230-2012-7710] 方程组的解为: 可看出,算得的该线性方程组的解向量为: [] 4 题目 (1)编制求第一型3次样条插值函数的通用程序; (2)已知汽车门曲线型值点的数据如下: i 0 1 2 3 4 5 6 7 8 9 10 Xi 0 1 2 3 4 5 6 7 8 9 10 Yi 端点条件为 ,用所编程序求车门的3次样条插值函数S(x),并打印出S(i+,i=0,1,…,9。 程序 clear digits(6); n=input('请输入节点数: n='); xn=zeros(1,n); yn=zeros(1,n); xn(1,: )=input('请输入节点坐标: '); yn(1,: )=input('请输入节点处函数值: '); dy0=input('请输入左边界条件: y’(x0)='); dyn=input('请输入右边界条件: y’(xn)='); %====================求d====================% d=zeros(n,1); h=zeros(1,n-1); f1=zeros(1,n-1); f2=zeros(1,n-2); fori=1: n-1 h(i)=xn(i+1)-xn(i); f1(i)=(yn(i+1)-yn(i))/h(i); end fori=2: n-1 f2(i)=(f1(i)-f1(i-1))/(xn(i+1)-xn(i-1)); d(i)=6*f2(i); end d(i)=6*(f1 (1)-dy0)/h (1); d(n)=6*(dyn-f1(n-1))/h(n-1); %====================求Mi====================% A=zeros(n); miu=zeros(1,n-2); lamda=zeros(1,n-2); fori=1: n-2 miu(i)=h(i)/(h(i)+h(i+1)); lamda(i)=1-miu(i); end A(1,2)=1; A(n,n-1)=1; fori=1: n A(i,i)=2; end fori=2: n-1 A(i,i-1)=miu(i-1); A(i,i+1)=lamda(i-1); end M=A\d; %====================回代求插值函数====================% symsx; fori=1: n-1; Sx(i)=collect(yn(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-xn(i))+M(i)/2*(x-xn(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-xn(i))^3); Sx(i)=vpa(Sx(i),6); end S=zeros(1,n-1); fori=1: n-1 x=xn(i)+; S(i)=yn(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-xn(i))+M(i)/2*(x-xn(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-xn(i))^3; end %====================打印结果====================% disp('S(x)='); fori=1: n-1 formatshort; fprintf('%s(%d disp('======================================================================'); end disp('S(i+') disp('ix(i+S(i+'); fori=1: n-1 fprintf('%d%.5f%.5f\n',i,xn(i)+,S(i)); end 运行结果 >>P195T37 请输入节点数: n=11 请输入节点坐标: [012345678910] 请输入节点处函数值: [] 请输入左边界条件: y’(x0)= 请输入右边界条件: y’(xn)= S(x)= *x+*x^2-*x^3+(0 ====================================================================== *x-*x^2-*x^3+(1 ====================================================================== *x-*x^2-*x^3+(2 ====================================================================== *x^2-*x-*x^3+(3 ====================================================================== *x-*x^2+*x^3-(4 ====================================================================== *x^2-*x-*x^3+(5 ====================================================================== *x-*x^2+*x^3-(6 ====================================================================== *x^2-*x-*x^3+(7 ====================================================================== *x-*x^2+*x^3-(8 ====================================================================== *x-*x^2+*x^3-(9 ====================================================================== S(i+ ix(i+S(i+ 1 2 3 4 5 6 7 8 9 10 5 题目 用Romberg求积法计算积分 的近似值,要求误差不超过 。 程序 %被积函数m文件: functionFx=fx(x) Fx=1/(1+100*x*x); end %Romberg求积法计算积分的通用程序 functionRomberg() clear; a=input('请输入积分下限: a='); b=input('请输入积分上限: b='); eps=input('请输入允许精度: eps='); %========计算Tn========% functionTn=T(n) Tn=0; h=(b-a)/n; x=zeros(1,n+1); fork=1: n+1 x(k)=a+(k-1)*h; end forj=1: n Tn=Tn+h*(fx(x(j))+fx(x(j+1)))/2; end end %========计算Sn========% functionSn=S(n) Sn=4/3*T(2*n)-1/3*T(n); end %========计算Cn========% functionCn=C(n) Cn=16/15*S(2*n)-1/15*S(n); end %========计算Rn========% functionRn=R(n) Rn=64/63*C(2*n)-1/63*C(n); end %========计算满足允许精度的Rn,并打印输出========% i=1; flag=1; whileflag==1 ifabs(R(2^i)-R(2^(i-1)))/255 flag=0; end i=i+1; end fprintf('该积分的值为: %f\n',R(2^(i-1))); end 运行结果 >>Romberg 请输入积分下限: a=-1 请输入积分上限: b=1 请输入允许精度: eps=*10^-7 该积分的值为: 结果分析 手动化简该定积分并最终求得的值为: ,误差限为: ,可见,程序完成了计算要求。 6 题目 常微分方程初值问题数值解 (1)编制RK4方法的通用程序; (2)编制AB4方法的通用程序(由RK4提供初值); (3)编制AB4-AM4预测校正方法通用程序(由RK4提供初值); (4)编制带改进的AB4-AM4预测校正方法通用程序(由RK4提供初值); (5)对于初值问题 取步长h=,应用 (1)-(4)中的四种方法进行计算,并将计算结果和精确解 作比较; (6)通过本上机题,你能得到哪些结论 程序 %f(x,y)函数m文件: functionFXY=fxy(x,y) FXY=-x*x*y*y; end %精确解y(x)函数m文件: functionFX=fx(x) FX=3/(1+x*x*x); end %RK4法通用程序 functionRK4() clear; x (1)=input('请输入初始x值: x0='); y (1)=input('请输入初值条件: y(x0)='); N=input('请输入计算步长: N='); h=input('请输入步长: h='); fori=1: N-1 x(i+1)=x(i)+h; k1=fxy(x(i),y(i)); k2=fxy(x(i)+*h,y(i)+*h*k1); k3=fxy(x(i)+*h,y(i)+*h*k2); k4=fxy(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4); end disp('ixiyiy(xi)y(xi)-yi'); disp('-------------------------------------------------------------'); fori=1: N fprintf('%d%f%f%f%f\n',i,x(i),y(i),fx(x(i)),fx(x(i))-y(i)); disp('-------------------------------------------------------------'); end end %AB4法通用程序 functionAB4() clear; x (1)=input('请输入初始x值: x0='); y (1)=input('请输入初值条件: y(x0)='); N=input('请输入计算步长: N='); h=input('请输入步长: h='); fori=1: N-1 x(i+1)=x(i)+h; k1=fxy(x(i),y(i)); k2=fxy(x(i)+*h,y(i)+*h*k1); k3=fxy(x(i)+*h,y(i)+*h*k2); k4=fxy(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4); end fori=4: N-1 y(i+1)=y(i)+h/24*(55*fxy(x(i),y(i))-59*fxy(x(i-1),y(i-1))+37*fxy(x(i-2),y(i-2))-9*fxy(x(i-3),y(i-3))); end disp('ixiyiy(xi)y(xi)-yi'); disp('-------------------------------------------------------------'); fori=1: N fprintf('%d%f%f%f%f\n',i,x(i),y(i),fx(x(i)),fx(x(i))-y(i)); disp('-------------------------------------------------------------'); end end %AB4-AM4预测校正法通用程序 functionAB4AM4() clear; x (1)=input('请输入初始x值: x0='); y (1)=input('请输入初值条件: y(x0)='); N=input('请输入计算步长: N='); h=input('请输入步长: h='); fori=1: N-1 x(i+1)=x(i)+h; k1=fxy(x(i),y(i)); k2=fxy(x(i)+*h,y(i)+*h*k1); k3=fxy(x(i)+*h,y(i)+*h*k2); k4=fxy(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4); end fori=4: N-1 yp(i+1)=y(i)+h/24*(55*fxy(x(i),y(i))-59*fxy(x(i-1),y(i-1))+37*fxy(x(i-2),y(i-2))-9*fxy(x(i-3),y(i-3))); y(i+1)=y(i)+h/24*(9*fxy(x(i+1),yp(i+1))+19*fxy(x(i),y(i))-5*fxy(x(i-1),y(i-1))+fxy(x(i-2),y(i-2))); end disp('ixiyiy(xi)y(xi)-yi'); disp('-------------------------------------------------------------'); fori=1: N fprintf('%d%f%f%f%f\n',i,x(i),y(i),fx(x(i)),fx(x(i))-y(i)); d
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 东南大学 数值 分析 上机 作业 MATLAB