东南大学数值分析编程作业

1、数值分析编程作业姓 名：卢一鸣学 号：151106 指导老师：吴宏伟目录P20.第17题：舍入误差与有效数3P56.第20题：Newton迭代法5P126.第39题：列主元Gauss消去法13P127.第40题：逐次超松弛迭代法16P195.第37题：3次样条插值函数21P256.第23题：重积分的计算23P307.第23题：常微分方程初值问题数值解25P346.第10题：抛物方程Crank-Nicolson格式30注：源程序放在对应的文件夹里，如第一题放在“P20.第17题：舍入误差与有效数”文件夹里；P20.第17题：舍入误差与有效数设（1）编制按从小到大的顺序，计算的通用程序；（2）编制

2、按从小到大的顺序，计算的通用程序；（3）按两种顺序分别计算，并指出有效位数（编制程序时用单精度）；（4）通过本上机题你明白了什么？答（1）（2）程序：程序放在computation.m, goon.m文件里，运行程序时只需输入命令main即可。%此函数为main函数%此函数第一个for循环N按从小到大顺序disp(blanks(0),本函数三种SN值：);disp(blanks(10),N按从小到大);disp(blanks(10),N按从大到小);disp(blanks(10),真实值);for N=100 10000 1000000 w=sprintf(%s%d,N=,N); disp(b

3、lanks(15),w); computation(N);endgoon();function a = computation( N )%本函数计算三种SN值：N按从小到大，N按从大到小，真实值%第二个for循环N按从大到小顺序%第二个for循环N按从大到小顺序SN=0;SN=single(SN);%N从小到大的顺序for i=2:N SN=SN+1/(i*i-1);enda=sprintf(%5.20f,SN);disp(blanks(10),a);SN=0;SN=single(SN);%N从大到小的顺序for i=N:-1:2 SN=SN+1/(i*i-1);enda=sprintf(%5

4、.20f,SN);disp(blanks(10),a);%真实值b=0.5*(1.5-1/N-1/(N+1);b=sprintf(%5.20f,SN);disp(blanks(10),b);endfunction a = goon( )%UNTITLED3 Summary of this function goes here% Detailed explanation goes herewhile 1disp(现在你可以尝试其他的N值,如想结束程序，请输入N=1);N=input(请输入N：);if N=1 return;else w=sprintf(%s%d,N=,N); disp(blan

5、ks(15),w); computation(N);endend（3）运行结果：N=100时，S1有6位有效数字，S2都是有效数字N=10000时，S1有4位有效数字，S2都是有效数字N=1000000时，S1有3位有效数字，S2都是有效数字（4）体会：当N从小到大变化时，SN的项从大到小，反之SN的项从小到大。从运算结果可见，第一种计算结果总是比第2、3种计算结果小，这是由于大数吃小数的原因（计算机表示一个数值是用字节表示的，当程序进行计算时，先将第一项和第二项相加，然后再加第三项所以加到后面会由于项的值很小，不能加到这个字节里，所以就被忽略了。）第2、3种计算结果一样，说明当N从大到小变化

6、时，可以避免大数吃小数，这在我们进行程序设计时，应尤其注意。P56.第20题：Newton迭代法（1）给定初值及容许误差，编制Newton法解方程根的通用程序。（2）给定方程 ，易知其有三个根，。由Newton方法的局部收敛性可知存在，当时Newton迭代序列收敛于根，试确定尽可能大的；试取若干初始值，观察当时Newton序列是否收敛以及收敛于哪一个根。（3）通过本上机题，你明白了什么？答：（1）通用程序在newton.m文件里。程序执行前在f.m和df.m文件中输入和，在命令行中输入main，并按要求即能获得结果。%此函数为main函数x0=input(请输入x0:);delta=5e-4;

7、newton(f,df,x0,delta ); %delta=search(f,df,delta )function y = f( x )%f(x)y=x3/3-x;endfunction x1,err,k,y = newton(f,df,x0,delta )%f是非线性函数%df是f的微商%x0是初始值%delta是给定允许误差%p1是牛顿法求得的方程的近似值%err是x0的误差估计%y=f(x1)%此程序没有求出f的微分，必须自己手动求出disp(blanks(3);disp(k表示迭代次数，xk表示根，err表示误差，y表示函数值)disp(k,blanks(8),xk,blanks(1

8、3),err,blanks(15),y);k=0;w=sprintf(%dt%5.10f,k,x0); disp(w);while 1 k=k+1; x1=x0-feval(f,x0)/feval(df,x0); err=abs(x1-x0); x0=x1; y=feval(f,x1); if(errdelta)|(y=0), w=sprintf(%dt%5.10ft%et%5.10f,k,x0,err,y); disp(w); break endendfunction delta = search( f,df,delta )%寻找尽可能大的delta flag=1;k=1;x0=0;whil

9、e flag=1 delta=k*5e-4; x0=delta; k=k+1; m=0; flag1=1; while flag1=1 & m=103 x1=x0-feval(f,x0)/feval(df,x0); if abs(x1-x0)=5e-4 flag=0; endendendfunction y = df( x )%f(x)y=x2-1; end（2）取容许误差为5e-4，通过search.m文件（执行命令delta=search(f,df,delta )）找到尽可能大的。 。区间上取-1000，-100，-50，-30，-10，-7，-5，-2，-1.5，-1.1，可见这些初始值

10、都收敛于。在区间即区间（-1，-0.7750）上取-0.99，-0.9，-0.85，-0.8，-0.7755， 可见有些初始值收敛于，而有些初始值收敛于。在区间即区间（-0.7750，0.7750）上，由search.m的运行过程表明，在整个区间上均收敛于。在区间即区间（0.7750，1）上取0.7755，0.8，0.85，0.9，0.99，可见有些初始值收敛于，而有些初始值收敛于。区间上取1000，100，50，30，10，7，5，2，1.5，1.1，可见初始值都收敛于。（3）综上所述：(-,-1)区间收敛于-1.73205,(-1,)区间一部分收敛于1.73205,一部分收敛于-1.732

11、05,(-,)区间收敛于0,(,1)区间类似于(-1,)区间，一部分收敛于1.73205,一部分收敛于-1.73205 ，(1,)收敛于1.73205。通过本上机题，明白了对于多根方程，Newton法求方程根时，迭代序列收敛于某一个根有一定的区间限制，在一个区间上，可能会局部收敛于不同的根。所以当我们求解这种方程时，最后先通过matlab的画图功能把曲线画出来，大概掌握零点的个数、位置以判断迭代出来的解是否正确和完整。总的来说，初始值距离收敛值越近，迭代次数越小，但是当差距很远（如1000）时，迭代次数也不是太多，大概20次，所以newton迭代法只需较少的迭代次数就能得到收敛值。P126.第

12、39题：列主元Gauss消去法对于某电路的分析，归结为求解线性方程组，其中（1）编制解n阶线性方程组的列主元Gauss消去法的通用程序；（2）用所编程序解线性方程组，并打印出解向量，保留5位有效数字；（3）本题编程之中，你提高了哪些编程能力？答：（1）通用程序在Gauss.m文件里。程序执行前在Aandb.m文件中输入系数矩阵A和右端向量，执行main命令，即能获得结果。%此函数为main函数clear allA,b=Aandb();x,flag =Gauss( A,b );function x,flag =Gauss( A,b )%求线性方程组的列主元Gauss消去法，其中，%A为方程组的系

13、数矩阵；%b为方程组的右端项；%x为方程组的解；%det为系数矩阵A的行列式的值；%flag为指标向量，flag=failure表示计算失败，flag=OK表示计算成功 A,bn,m=size(A);nb=length(b);%当方程组行与列的维数不相等时，停止计算，并输出出错信息 if n=m error(The rows and columns of matrix A must be equal!); return; end %当方程组与右端项的维数不匹配时，停止计算，并输出出错信息 if m=nb error(The columns of A must be equal the leng

14、th of b!) return; end % 开始计算，先赋初值 flag=OK;x=zeros(n,1); for k=1:n-1 %选主元 max1=0; for i=k:n if abs(A(i,k)max1 max1=abs(A(i,k);r=i; end end if max1k for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z; end %肖元过程 for i=k+1:n l=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-l*A(k,j); end b(

15、i)=b(i)-l*b(k); end end%回代过程 if abs(A(n,n)1e-10 flag=failure;return; end for k=n:-1:1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); end disp(x=); fprintf(%.5gn,x) endfunction A,b = Aandb( )%系数矩阵和右端向量% Detailed explanation goes hereA=31 -13 0 0 0 -10 0 0 0; -13 35 -9 0 -11 0 0 0 0; 0 -9 31

16、 -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;%A=2 -4 6;% 4 -9 2;% 1 -1 3;%b=3;5;4; end（2）首先将R和V的内容输入到Aandb.m文件中。输入命令main即得结果。（3）本题编程过程中，由于是一段通用程序，要考虑到每一种可能，比如有计算成

17、功的，有计算不成功的。矩阵数据较多，在输出界面上最好能够再一次显示矩阵内容，方便用户查看矩阵输入是否正确。要善于利用循环，掌握好每一个变量的意义，顺序，一个小地方出错则很可能导致整个程序运行不起来。所以说编程是一项细活。P127.第40题：逐次超松弛迭代法（1）编制解n阶线性方程组的SOR方法的通用程序（要求）；（2）对于第39题中所给的线性方程组，取松弛因子，容许误差，打印松弛因子、迭代次数、最佳松弛因子及解向量。答：（1）通用程序在SOR.m文件里。程序执行前在Aandb.m文件中输入系数矩阵A和右端向量，执行main命令，即能获得结果。%此函数为main函数clear allA,b=Aa

18、ndb();delta=0.5e-5;max=10000; %如果迭代次数超过10000，失败disp( w );%x,k,flag =SOR( A,b,delta,1.96,max);for i=1:99 w=i/50; x,k,flag =SOR( A,b,delta,w,max); a=sprintf(%1.2ft%1.0f,w,k); disp(a); q(i)=w; c(i)=k;endm,n=min(c);x,k,flag =SOR( A,b,delta,q(n),max);a=sprintf(%s%1.2fn%s, 最佳松弛因子：,q(n),解向量);disp(a);fprint

19、f(%5.10gn,x);function x,k,flag = SOR( A,b,delta,w,max )%求解线性方程组的迭代法，其中%A为方程组的系数矩阵%b为方程组的右端项%delta为精度要求，缺省值为1e-5%max1为最大迭代次数，缺省值100%w为超松弛因子，缺省值为1%x为方程组的解%k为迭代次数%flag为指标变量， flag=OK！表示迭代收敛到指标要求% flag=failure！表示迭代失败 if nargin5 max=10000;endif nargin4 w=1;endif nargin3 delta=1e-5;endn=length(A);k=0;x=zer

20、os(n,1);y=zeros(n,1);flag=OK!;while 1 y=x; for i=1:n z=b(i); for j=1:n if j=i z=z-A(i,j)*x(j); end end if abs(A(i,i)1e-10|k=max flag=fail; return; end z=z/A(i,i);x(i)=(1-w)*x(i)+w*z; end if norm(y-x,inf)1 T(k-1,2)=4/3*T(k,1)-1/3*T(k-1,1); end if k2 T(k-2,3)=16/15*T(k-1,2)-1/15*T(k-2,2); end if k3 T(

21、k-3,4)=64/63*T(k-2,3)-1/63*T(k-3,3); end if k4 err=1/255*abs(T(k-3,4)-T(k-4,4); if(errdelta)|k=20 w=sprintf(%s%2.8f, 满足精确度条件的近似值为：,T(k-3,4); disp(w); break end end endm,n=size(T);w=sprintf(%st%st%st%st%st%s,k,2k, T(2k), S(2k), C(2k), R(2k);disp(w);for i=1:mfprintf(%dt%dt%2.8ft%2.8ft%2.8ft%2.8fn,i,2i

22、,T(i,1),T(i,2),T(i,3),T(i,4) endfunction y = f( x,y )%f(x)y=tan(x2+y2);end function T = complexinteg(a,b,m,c,d,n )%本函数用复化梯形公式求重积分%f是被积函数,%a,b分别为x的上下限,c,d为y的上下限，%a，b,c，c分别划分成m，n份 %a=0;b=3;m=3;d=6;c=4;n=2; h=(b-a)/m; x=a:h:b; k=(c-d)/n; y=d:k:c;%x,y=meshgrid(x,y);T=0;%角点c=feval(f,x(1),y(1);T=T+feval(f

23、,x(1),y(1)+feval(f,x(1),y(n+1)+feval(f,x(m+1),y(1)+feval(f,x(m+1),y(n+1);%边点for i=1 m+1 for j=2:n T=T+2*feval(f,x(i),y(j); endendfor j=1 n+1 for i=2:m T=T+2*feval(f,x(i),y(j); endend %内部点for i=2:m for j=2:n T=T+4*feval(f,x(i),y(j); endendT=-T*h*k/4; （2）计算结果如下。说明：显示中有许多0，这是因为我用的矩阵记录每个值，没有算的位置程序就默认为0，

24、这忽略。 体会：此题程序不算复杂，只要有一些基本功，再加上掌握算法，就可以编写出来，这难道是之前几道题的认真编程产生的后果。学习永无止境。P307.第23题：常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法通用程度（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；（5）对于初值问题取步长h=0.1，应用（1）（4）中的四种方法进行计算，并将计算结果和精确解作比较。（6）通过本上机题，你能得到哪些结论？答：RK4：AB4：AB4-AM4：改进的AB4-AM

25、4：（1）RK4，AB4，AB4-AM4，改进的AB4-AM4算法放在文件RK4.m,AB4.m, AB4-AM4.m, modiAB4-AM4.m,accuracy.m文件里，根据已给解析解方程给出离散点处的y值，在文件accuracy.m。程序执行输入main命令，即能获得结果。说明：在程序执行前， f.m里存放f(x,y)；f1.m里存放解析解的表达式。%此函数为main函数rk4=RK4(f,0,1.5,3,15);ab4=AB4( f,0,1.5,rk4(1,2); %RK4提供初始值rk4(2,2),rk4(3,2),rk4(4,2),15 );ab4am4=AB4AM4( f,0

26、,1.5,rk4(1,2); rk4(2,2),rk4(3,2),rk4(4,2),15 );modiab4am4=modiAB4AM4( f,0,1.5,rk4(1,2); rk4(2,2),rk4(3,2),rk4(4,2),15 );acc=accuracy( f,0,1.5,15 ); %真实解format long;a=acc,rk4(:,2),ab4(:,2),ab4am4(:,2),modiab4am4(:,2);b=abs(rk4(:,2)-acc(:,2),abs(ab4(:,2)-acc(:,2),abs(ab4am4(:,2)-acc(:,2),abs(modiab4am

27、4(:,2)-acc(:,2); %误差显示fprintf(tt%sttttt%sttttt%stttt%sttttt%sttt%sn,xk, 精确解,RK4,AB4,AB4-AM4,改进的AB4-AM4);disp(a);fprintf(tt%sttt%stt%stt%sn, RK4的误差,AB4的误差,AB4-AM4的误差,改进的AB4-AM4的误差);disp(b);function y = f( x,y )%f(x)y=-x2*y2;endfunction y = f1(x )%真实解y=3/(1+x3);endfunction R = RK4( f,a,b,y0,N )%RK4算法%

28、y=f（x，y）%a,b左右端点%N为迭代步长%h为步长%y0为初值h=(b-a)/N;X=zeros(1,N+1);Y=zeros(1,N+1);X=a:h:b;Y(1)=y0;for j=1:N k1=h*feval(f,X(j),Y(j); k2=h*feval(f,X(j)+h/2,Y(j)+k1/2); k3=h*feval(f,X(j)+h/2,Y(j)+k2/2); k4=h*feval(f,X(j)+h,Y(j)+k3); Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endR=X Y;function R = AB4( f,a,b,y0,y1,y2,y3,N

29、 )%AB4算法%y=f（x，y）%a,b左右端点%N为迭代步长%h为步长%y0,y1,y2,y3为初值,RK4提供h=(b-a)/N;X=zeros(1,N+1);Y=zeros(1,N+1);X=a:h:b;Y(1)=y0;Y(2)=y1;Y(3)=y2;Y(4)=y3;for j=4:N Y(j+1)=Y(j)+h/24*(55*feval(f,X(j),Y(j)-59*feval(f,X(j-1),Y(j-1)+37*feval(f,X(j-2),Y(j-2)-9*feval(f,X(j-3),Y(j-3);endR=X Y;endfunction R = AB4AM4( f,a,b,

30、y0,y1,y2,y3,N )%AB4-AM4算法%y=f（x，y）%a,b左右端点%N为迭代步长%h为步长%y0,y1,y2,y3为初值,RK4提供h=(b-a)/N;X=zeros(1,N+1);Y=zeros(1,N+1);X=a:h:b;Y(1)=y0;Y(2)=y1;Y(3)=y2;Y(4)=y3;for j=4:N Y(j+1)=Y(j)+h/24*(55*feval(f,X(j),Y(j)-59*feval(f,X(j-1),Y(j-1)+37*feval(f,X(j-2),Y(j-2)-9*feval(f,X(j-3),Y(j-3); Y(j+1)=Y(j)+h/24*(9*feval(f,X(j+1),Y(j+1)+19*feval(f,X(j),Y(j)-5*feval(f,X(j-1),Y(j-1)+feval(f,X(j-2),Y(j-2);endR=X Y;end function R = modiAB4AM4(