东南大学数值分析上机题作业matlab版(共18页).doc

《东南大学数值分析上机题作业matlab版(共18页).doc》由会员分享，可在线阅读，更多相关《东南大学数值分析上机题作业matlab版(共18页).doc（18页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上2015.上机作业题报告 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);for i=2:N Snl2s=Snl2s+1/(i*i-1);endfor i=N:-1:2 Sns2l=Sns2l+1/(i*i-1);endfprintf(精

2、确值为： %fn,Ac);fprintf(从大到小的顺序累加得SN=%fn,Snl2s);fprintf(从小到大的顺序累加得SN=%fn,Sns2l);disp(=);程序 P20T17请输入N值:102精确值为： 从大到小的顺序累加得SN=从小到大的顺序累加得SN= P20T17请输入N值:104精确值为： 从大到小的顺序累加得SN=运行结果从小到大的顺序累加得SN= P20T17请输入N值:106精确值为： 从大到小的顺序累加得SN=从小到大的顺序累加得SN=结果分析按从大到小的顺序，有效位数分别为：6，4，3。按从小到大的顺序，有效位数分别为：5，6，6。可以看出，不同的算法造成的误差

3、限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2题目（1）给定初值及容许误差，编制牛顿法解方程f(x)=0的通用程序。（2）给定方程,易知其有三个根由牛顿方法的局部收敛性可知存在当时，Newton迭代序列收敛于根x2*。试确定尽可能大的。试取若干初始值，观察当时Newton序列的收敛性以及收敛于哪一个根。（3）通过本上机题，你明白了什么？程序函数m文件：function Fu=fu(x)Fu=x3/3-x;end函数m文件：function

4、 Fu=dfu(x)Fu=x2-1;end用Newton法求根的通用程序clear;x0=input(请输入初值x0:);ep=input(请输入容许误差:);flag=1;while flag=1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)epflag=0; endx0=x1;endfprintf(方程的一个近似解为：%fn,x0);寻找最大值的程序：cleareps=input(请输入搜索精度：);ep=input(请输入容许误差：);flag=1;k=0;x0=0;while flag=1 sigma=k*eps;x0=sigma; k=k+1; m=0; f

5、lag1=1; while flag1=1 & m=103 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0; endendfprintf(最大的sigma值为：%fn,sigma);运行结果（1）寻找最大的值。算法为：将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。运行，得到在不同的搜索精度下的最大sigma值。 Find请输入搜索精度：10-6请输入容许误差：10-6最大的sigma值为： Find请输入搜索精度：10-4请输入容许误差：10-6最大的sigma

6、值为： Find请输入搜索精度：10-2请输入容许误差：10-6最大的sigma值为：（2）运行在内取初值，运行结果如下：X0Xk-1000-500-100-10-5可见，在区间内取初值，Newton序列收敛，且收敛于根。在内取初值，运行结果如下：X0Xk可见，在内取初值，Newton序列收敛，且收敛于根。在内内取初值，运行结果如下：X0Xk可见，在内取初值，Newton序列收敛，且收敛于根0。在内取初值，运行结果如下：X0Xk可见，在内取初值，Newton序列收敛，且收敛于根在内取初值，运行结果如下：X0Xk5101005001000可见，在内取初值，Newton序列收敛，且收敛于根 3题目

7、对于某电路的分析，归结为求解线性方程组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;for i=1:n-1 maximum,index=max(abs(C(i:n,i); index=index+i-1; T=C(index,:); C(index,:)=C(i,:)

8、; C(i,:)=T; for k=i+1:n if C(k,i)=0 C(k,:)=C(k,:)-C(k,i)/C(i,i)*C(i,:); end endend%回代求解x=zeros(n,1);x(n)=C(n,n+1)/C(n,n);for i=n-1:-1:1 x(i)=(C(i,n+1)-C(i,i+1:n)*x(i+1:n,1)/C(i,i);enddisp(方程组的解为:);fprintf(%.5gn,x);运行结果运行程序，输入系数矩阵和方程组右端列向量。运行过程与结果如下图所示： P126T39请输入线性方程组阶数: n=4请输入系数矩阵：A= 0 0; 0;0 ;0 0

9、请输入线性方程组右端向量：b= 方程组的解为:2495 P126T39请输入线性方程组阶数: n=9请输入系数矩阵：A=31 -13 0 0 0 -10 0 0 0;-13 35 -9 0 -11 0 0 0 0;0 -9 31 -10 0 0 0 0 0;0 0 -10 79 -30 0 0 0 -9;0 0 0 -30 57 -7 0 -5 0;0 0 0 0 -7 47 -30 0 0;0 0 0 0 0 -30 41 0 0;0 0 0 0 -5 0 0 27 -2;0 0 0 -9 0 0 0 -2 29请输入线性方程组右端向量：b=-15 27 -23 0 -20 12 -7 7

10、10方程组的解为:可看出，算得的该线性方程组的解向量为： 4题目（1）编制求第一型3次样条插值函数的通用程序；（2）已知汽车门曲线型值点的数据如下：i012345678910Xi012345678910Yi端点条件为，用所编程序求车门的3次样条插值函数S(x),并打印出S(i+，i=0,1，9。程序cleardigits(6);n=input(请输入节点数：n=);xn=zeros(1,n);yn=zeros(1,n);xn(1,:)=input(请输入节点坐标：);yn(1,:)=input(请输入节点处函数值：);dy0=input(请输入左边界条件：y（x0）=);dyn=input(请

11、输入右边界条件：y（xn）=);%=求d=%d=zeros(n,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n-1 h(i)=xn(i+1)-xn(i); f1(i)=(yn(i+1)-yn(i)/h(i);endfor i=2:n-1 f2(i)=(f1(i)-f1(i-1)/(xn(i+1)-xn(i-1); d(i)=6*f2(i);endd(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);lamd

12、a=zeros(1,n-2);for i=1:n-2 miu(i)=h(i)/(h(i)+h(i+1); lamda(i)=1-miu(i);endA(1,2)=1;A(n,n-1)=1;for i=1:n A(i,i)=2;endfor i=2:n-1 A(i,i-1)=miu(i-1); A(i,i+1)=lamda(i-1);endM=Ad;%=回代求插值函数=%syms x;for i=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(

13、i)*(x-xn(i)3); Sx(i)=vpa(Sx(i),6);endS=zeros(1,n-1);for i=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)=);for i=1:n-1 format short; fprintf( %s (%dx P195T37请输入节点数：n=11请输入节点坐标：0 1 2 3 4 5 6 7 8 9 10请输入节点处函数值： 请输入

14、左边界条件：y（x0）=请输入右边界条件：y（xn）=S(x)= *x + *x2 - *x3 + (0x1)= *x - *x2 - *x3 + (1x2)= *x - *x2 - *x3 + (2x3)= *x2 - *x - *x3 + (3x4)= *x - *x2 + *x3 - (4x5)= *x2 - *x - *x3 + (5x6)= *x - *x2 + *x3 - (6x7)= *x2 - *x - *x3 + (7x8)= *x - *x2 + *x3 - (8x9)= *x - *x2 + *x3 - (9x10)=S(i+ i x(i+ S(i+ 1 2 3 4 5

15、6 7 8 9 10 5题目用Romberg求积法计算积分的近似值，要求误差不超过。程序%被积函数m文件：function Fx=fx(x)Fx=1/(1+100*x*x);end%Romberg求积法计算积分的通用程序function Romberg()clear;a=input(请输入积分下限：a=);b=input(请输入积分上限：b=);eps=input(请输入允许精度：eps=);%=计算Tn=%function Tn=T(n)Tn=0;h=(b-a)/n;x=zeros(1,n+1);for k=1:n+1 x(k)=a+(k-1)*h;endfor j=1:n Tn=Tn+h*

16、(fx(x(j)+fx(x(j+1)/2;endend%=计算Sn=%function Sn=S(n)Sn=4/3*T(2*n)-1/3*T(n);end%=计算Cn=%function Cn=C(n)Cn=16/15*S(2*n)-1/15*S(n);end%=计算Rn=%function Rn=R(n)Rn=64/63*C(2*n)-1/63*C(n);end%=计算满足允许精度的Rn，并打印输出=%i=1;flag=1;while flag=1if abs(R(2i)-R(2(i-1)/255 Romberg请输入积分下限：a=-1请输入积分上限：b=1请输入允许精度：eps=*10-7

17、该积分的值为：结果分析手动化简该定积分并最终求得的值为：，误差限为：，可见，程序完成了计算要求。 6题目常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；（5）对于初值问题取步长h=，应用（1）-（4）中的四种方法进行计算，并将计算结果和精确解作比较；（6）通过本上机题，你能得到哪些结论？程序%f（x，y）函数m文件：function FXY=fxy(x,y)FXY=-x*x*y*y;end%精确解

18、y（x）函数m文件：function FX=fx(x)FX=3/(1+x*x*x);end%RK4法通用程序function RK4()clear;x(1)=input(请输入初始x值：x0=);y(1)=input(请输入初值条件：y（x0）=);N=input(请输入计算步长：N=);h=input(请输入步长：h=);for i=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

19、)+h/6*(k1+2*k2+2*k3+k4);end disp(i xi yi y(xi) y(xi)-yi); disp(-);for i=1:N fprintf(%d %f %f %f %fn,i,x(i),y(i),fx(x(i),fx(x(i)-y(i); disp(-);endend%AB4法通用程序function AB4()clear;x(1)=input(请输入初始x值：x0=);y(1)=input(请输入初值条件：y（x0）=);N=input(请输入计算步长：N=);h=input(请输入步长：h=);for i=1:N-1 x(i+1)=x(i)+h; k1=fxy(

20、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);endfor i=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(i xi yi y(xi) y(xi)-yi); disp(-);for i=1:N fprintf(%

21、d %f %f %f %fn,i,x(i),y(i),fx(x(i),fx(x(i)-y(i); disp(-);endend%AB4-AM4预测校正法通用程序function AB4AM4()clear;x(1)=input(请输入初始x值：x0=);y(1)=input(请输入初值条件：y（x0）=);N=input(请输入计算步长：N=);h=input(请输入步长：h=);for i=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(

22、x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);endfor i=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(i xi yi y(xi) y(xi)-yi); disp(

23、-);for i=1:N fprintf(%d %f %f %f %fn,i,x(i),y(i),fx(x(i),fx(x(i)-y(i); disp(-);endend%带改进的AB4-AM4预测校正方法function AB4AM4plus()clear;x(1)=input(请输入初始x值：x0=);y(1)=input(请输入初值条件：y（x0）=);N=input(请输入计算步长：N=);h=input(请输入步长：h=);for i=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

24、(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);endfor i=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); yc(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); y(i+1)=251

25、/270*yc(i+1)+19/270*yp(i+1);end disp(i xi yi y(xi) y(xi)-yi); disp(-);for i=1:N fprintf(%d %f %f %f %fn,i,x(i),y(i),fx(x(i),fx(x(i)-y(i); disp(-);endend运行结果（1）RK4法 RK4请输入初始x值：x0=0请输入初值条件：y（x0）=3请输入计算步数：N=15请输入步长：h=i xi yi y(xi) y(xi)-yi-1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 -13 -14 -15 -（2）AB4法 AB4请输入初始x值：x0=0请输入初值条件：y（x0）=3请输入计算步数：N=15请输入步长：h=i xi yi y(xi) y(xi)-yi-1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 -13 -14 -15 -（3）AB4-AM4预测校