《东南大学数值分析编程作业》由会员分享，可在线阅读，更多相关《东南大学数值分析编程作业（42页珍藏版）》请在金锄头文库上搜索。
1、1数值分析编程作业姓 名：卢一鸣学 号：151106 2指导老师：吴宏伟3目录P20.第 17 题：舍入误差与有效数 .3P56.第 20 题：Newton 迭代法 .5P126.第 39 题：列主元 Gauss 消去法 .13P127.第 40 题：逐次超松弛迭代法 .16P195.第 37 题：3 次样条插值函数 .21P256.第 23 题：重积分的计算 .23P307.第 23 题：常微分方程初值问题数值解 .25P346.第 10 题：抛物方程 Crank-Nicolson 格式 .304注：源程序放在对应的文件夹里，如第一题放在“P20.第 17 题：舍入误差与有效数”文件夹里；P
2、20.第 17 题：舍入误差与有效数设 2113,-2N+NjS其 精 确 值 为 （ ） 。（1 ）编制按从小到大的顺序 ，计算 的通用程序；221-13-SNS（2 ）编制按从小到大的顺序 ，计算 的通用程序；222+N（ ）（3 ）按两种顺序分别计算 ， ， ，并指出有效位数（编制程序时用单精度） ；210S4610（4 ）通过本上机题你明白了什么？答（1 ） （2）程序：程序放在 computation.m, goon.m 文件里，运行程序时只需输入命令 main即可。%此函数为main函数%此函数第一个for循环N按从小到大顺序disp(blanks(0),本函数三种SN值：);di
3、sp(blanks(10),N按从小到大);disp(blanks(10),N按从大到小);disp(blanks(10),真实值 );for N=100 10000 1000000w=sprintf(%s%d,N=,N);disp(blanks(15),w);5computation(N);endgoon();function a = computation( N )%本函数计算三种SN值：N按从小到大，N按从大到小，真实值%第二个for 循环N按从大到小顺序%第二个for 循环N按从大到小顺序SN=0;SN=single(SN);%N从小到大的顺序for i=2:NSN=SN+1/(i*i
4、-1);enda=sprintf(%5.20f,SN);disp(blanks(10),a);SN=0;SN=single(SN);%N从大到小的顺序for i=N:-1:2SN=SN+1/(i*i-1);enda=sprintf(%5.20f,SN);disp(blanks(10),a);function a = goon( )%UNTITLED3 Summary of this function goes here% Detailed explanation goes herewhile 1disp(现在你可以尝试其他的N值, 如想结束程序，请输入N=5e-4 flag=0; functi
5、on y = df( x )%f(x)y=x2-1;end10endendend（2）取容许误差为 5e-4，通过 search.m 文件（执行命令 delta=search(f,df,delta )）找到尽可能大的 。0.75 。0(,1)(,.),(750,),(.750,1)(x区间 上取-1000，-100，-50 ，-30，-10，-7 ，-5，-2，-1.5 ，-1.1，可见这些初始值,都收敛于 。13x1112在区间 即区间（-1， -0.7750）上取-0.99 ，-0.9，-0.85，-0.8，-0.7755， 可见(1,)有些初始值收敛于 ，而有些初始值收敛于 。3x3x在
6、区间 即区间（-0.7750，0.7750）上，由 search.m 的运行过程表明，在整个(,)区间上均收敛于 。20x在区间 即区间（0.7750，1）上取 0.7755，0.8，0.85，0.9 ，0.99，可见有些初始(,)13值收敛于 ，而有些初始值收敛于 。13x3x区间 上取 1000，100，50，30，10 ，7 ，5，2 ，1.5，1.1，可见初始值都收敛于(1,)。3x141516（3 ） 综上所述：(-,-1)区间收敛于-1.73205,(-1,)区间一部分收敛于 1.73205,一部分收敛于-1.73205,(- ,)区间收敛于 0,(,1)区间类似于(-1,)区间，
7、一部分收敛于 1.73205,一部分收敛于-1.73205 ，(1,) 收敛于 1.73205。通过本上机题，明白了对于多根方程，Newton 法求方程根时，迭代序列收敛于某一个根有一定的区间限制，在一个区间上，可能会局部收敛于不同的根。所以当我们求解这种方程时，最后先通过 matlab 的画图功能把曲线画出来，大概掌握零点的个数、位置以判断迭代出来的解是否正确和完整。总的来说，初始值距离收敛值越近，迭代次数越小，但是当差距很远（如 1000）时，迭代次数也不是太多，大概 20 次，所以 newton 迭代法只需较少的迭代次数就能得到收敛值。P126.第 39 题：列主元 Gauss 消去法对
8、于某电路的分析，归结为求解线性方程组 ，其中RIV310100597390057054120909R (15,273,17,0)T TV（1 ）编制解 n 阶线性方程组 的列主元 Gauss 消去法的通用程序；Axb（2 ）用所编程序解线性方程组 ，并打印出解向量，保留 5 位有效数字；RI（3 ）本题编程之中，你提高了哪些编程能力？答：（1 ） 通用程序在 Gauss.m 文件里。程序执行前在 Aandb.m 文件中输入系数矩阵 A 和右端17向量，执行 main 命令，即能获得结果。%此函数为main函数clear allA,b=Aandb();x,flag =Gauss( A,b );f
9、unction x,flag =Gauss( A,b )%求线性方程组的列主元Gauss消去法，其中，%A为方程组的系数矩阵；%b为方程组的右端项；%x为方程组的解；%det为系数矩阵A的行列式的值；%flag为指标向量，flag= failure表示计算失败，flag=OK 表示计算成功A,bn,m=size(A);nb=length(b);%当方程组行与列的维数不相等时，停止计算，并输出出错信息if n=merror(The rows and columns of matrix A must be equal!);return;end %当方程组与右端项的维数不匹配时，停止计算，并输出出错
10、信息if m=nbfunction 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 -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;-
11、7;7;10;%A=2 -4 6;% 4 -9 2;% 1 -1 3;%b=3;5;4;18error(The columns of A must be equal the length of b!)return;end % 开始计算，先赋初值 flag=OK;x=zeros(n,1);for k=1:n-1%选主元max1=0;for i=k:nif abs(A(i,k)max1max1=abs(A(i,k);r=i;endendif max1kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;endend
12、19%肖元过程for i=k+1:nl=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-l*A(k,j);endb(i)=b(i)-l*b(k);endend%回代过程if abs(A(n,n)1function y = f( x,y )%f(x)y=tan(x2+y2);end33T(k-1,2)=4/3*T(k,1)-1/3*T(k-1,1);endif k2T(k-2,3)=16/15*T(k-1,2)-1/15*T(k-2,2);endif k3T(k-3,4)=64/63*T(k-2,3)-1/63*T(k-3,3);endif k4err=1/255*
13、abs(T(k-3,4)-T(k-4,4); if(errdelta)|k=20w=sprintf(%s%2.8f, 满足精确度条件的近似值为：,T(k-3,4);disp(w);breakendendendm,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,T(i,1),T(i,2),T(i,3),T(i,4) 34endfunction T = complexinteg(a,b
14、,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,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+1for j=2:nT=T+2*feval(f,x(i),y(j); endend35for j=1 n+1for i=2:mT=T+2*feval(f,x(i),y(j); endend %内部点for i=2:mfor j=2:nT=T+4*feval(f,x(i),y(j);endend
