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

1数值分析编程作业姓名：李红兵学号：130891指导老师：曹婉容2目录P20.第17题：舍入误差与有效数..........................................3P56.第20题：Newton迭代法................................................5P126.第39题：列主元Gauss消去法..................................13P127.第40题：逐次超松弛迭代法......................................16P195.第37题：3次样条插值函数.......................................22P256.第23题：重积分的计算..............................................24P307.第23题：常微分方程初值问题数值解......................26P346.第10题：抛物方程Crank-Nicolson格式...................313注：源程序放在对应的文件夹里，如第一题放在“P20.第17题：舍入误差与有效数”文件夹里；P20.第17题：舍入误差与有效数设2211311,--122NN+1NNjSj其精确值为（）。（1）编制按从小到大的顺序222111+++2-13-1N-1NS，计算NS的通用程序；（2）编制按从小到大的顺序222111+++N-1N-1-12-1NS（），计算NS的通用程序；（3）按两种顺序分别计算210S，410S，610S，并指出有效位数（编制程序时用单精度）；（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),'真实值']);forN=[100100001000000]w=sprintf('%s%d','N=',N);disp([blanks(15),w]);computation(N);endgoon();function[a]=computation(N)%本函数计算三种SN值：N按从小到大，N按从大到小，真实值%第二个for循环N按从大到小顺序%第二个for循环N按从大到小顺序SN=0;SN=single(SN);%N从小到大的顺序fori=2:NSN=SN+1/(i*i-1);enda=sprintf('%5.20f',SN);disp([blanks(10),a]);SN=0;function[a]=goon()%UNTITLED3Summaryofthisfunctiongoeshere%Detailedexplanationgoesherewhile1disp('现在你可以尝试其他的N值,如想结束程序，请输入N=1');N=input('请输入N：');ifN=1return;elsew=sprintf('%s%d','N=',N);disp([blanks(15),w]);computation(N);4SN=single(SN);%N从大到小的顺序fori=N:-1:2SN=SN+1/(i*i-1);enda=sprintf('%5.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]);endendend（3）运行结果：N=100时，S1有6位有效数字，S2都是有效数字N=10000时，S1有4位有效数字，S2都是有效数字N=1000000时，S1有3位有效数字，S2都是有效数字（4）体会：当N从小到大变化时，SN的项从大到小，反之SN的项从小到大。从运算结果可见，第一种计算结果总是比第2、3种计算结果小，这是由于大数吃小数的原因（计算机表示一个数值是用字节表示的，当程序进行计算时，先将第一项和第二项相加，然后再加第三项……所以加到后面会由于项的值很小，不能加到这个字节里，所以就被忽略5了。）第2、3种计算结果一样，说明当N从大到小变化时，可以避免大数吃小数，这在我们进行程序设计时，应尤其注意。P56.第20题：Newton迭代法（1）给定初值0x及容许误差，编制Newton法解方程()0fx根的通用程序。（2）给定方程3()/30fxxx，易知其有三个根13x，20x，33x。①由Newton方法的局部收敛性可知存在0，当0(,)x时Newton迭代序列收敛于根2x，试确定尽可能大的；②试取若干初始值，观察当0(,1),(1,),(,),(,1),(1,)x时Newton序列是否收敛以及收敛于哪一个根。（3）通过本上机题，你明白了什么？答：（1）通用程序在newton.m文件里。程序执行前在f.m和df.m文件中输入()fx和'()fx，在命令行中输入main，并按要求即能获得结果。%此函数为main函数x0=input('请输入x0:');delta=5e-4;newton('f','df',x0,delta);%delta=search('f','df',delta)function[y]=f(x)%f(x)y=x^3/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(13),'err',blanks(15),'y']);k=0;w=sprintf('%d\t%5.10f',k,x0);disp(w);while1k=k+1;6x1=x0-feval('f',x0)/feval('df',x0);err=abs(x1-x0);x0=x1;y=feval('f',x1);if(errdelta)|(y==0),w=sprintf('%d\t%5.10f\t%e\t%5.10f',k,x0,err,y);disp(w);breakendendfunction[delta]=search(f,df,delta)%寻找尽可能大的deltaflag=1;k=1;x0=0;whileflag==1delta=k*5e-4;x0=delta;k=k+1;m=0;flag1=1;whileflag1==1&&m=10^3x1=x0-feval('f',x0)/feval('df',x0);ifabs(x1-x0)5e-4flag1=0;endm=m+1;x0=x1;endifflag1==1||abs(x0)=5e-4flag=0;endendendfunction[y]=df(x)%f'(x)y=x^2-1;end（2）①取容许误差为5e-4，通过search.m文件（执行命令delta=search('f','df',delta)）找到尽可能大的0.7750。②0(,1),(1,0.7750),(7750,7750),(0.7750,1),(1,)x。区间(,1)上取-1000，-100，-50，-30，-10，-7，-5，-2，-1.5，-1.1，可见这些初始值都收敛于13x。78在区间(1,)即区间（-1，-0.7750）上取-0.99，-0.9，-0.85，-0.8，-0.7755，可见有些初始值收敛于13x，而有些初始值收敛于33x。9在区间(,)即区间（-0.7750，0.7750）上，由search.m的运行过程表明，在整个区间上均收敛于20x。在区间(,1)即区间（0.7750，1）上取0.7755，0.8，0.85，0.9，0.99，可见有些初始值收敛于13x，而有些初始值收敛于33x。10区间(1,)上取1000，100，50，30，10，7，5，2，1.5，1.1，可见初始值都收敛于33x。1112（3）综上所述：(-∞,-1)区间收敛于-1.73205,(-1,δ)区间一部分收敛于1.73205,一部分收敛于-1.73205,(-δ,δ)区间收敛于0,(δ,1)区间类似于(-1,δ)区间，一部分收敛于1.73205,一部分收敛于-1.73205，(1,∞)收敛于1.73205。13通过本上机题，明白了对于多根方程，Newton法求方程根时，迭代序列收敛于某一个根有一定的区间限制，在一个区间上，可能会局部收敛于不同的根。所以当我们求解这种方程时，最后先通过matlab的画图功能把曲线画出来，大概掌握零点的个数、位置以判断迭代出来的解是否正确和完整。总的来说，初始值距离收敛值越近，迭代次数越小，但是当差距很远（如1000）时，迭代次数也不是太多，大概20次，所以newton迭代法只需较少的迭代次数就能得到收敛值。P126.第39题：列主元Gauss消去法对于某电路的分析，归结为求解线性方程组RIV，其中3113000100001335901100000931100000000107930000900030577050000074730000000030410000005002720009000229R(15,27,23,0,20,12,7,7,10)TTV（1）编制解n阶线性方程组Axb的列主元Gauss消去法的通用程序；（2）用所编程序解线性方程组RIV，并打印出解向量，保留5位有效数字；（3）本题编程之中，你提高了哪些编程能力？答：（1）通用程序在Gauss.m文件里。程序执行前在Aandb.m文件中输入系数矩阵A和右端向量，执行main命令，即能获得结果。%此函数为main函数clearall[A,b]=Aandb();[x,flag]=Gauss(A,b);function[x,flag]=Gauss(A,b)%求线性方程组的列主元Gauss消去法，其中，%A为方程组的系数矩阵；%b为方程组的右端项；%x为方程组的解；%det为系数矩阵A的行列式的值；%flag为指标向量，flag=‘failure’表示计算失败，flag=‘OK’表示计算成功function[A,b]=Aandb()%系数矩阵和右端向量%DetailedexplanationgoeshereA=[31-13000-10000;-1335-90-110000;0-931-1000000;00-1079-30000-9;000-3057-70-50;14A,b[n,m]=size(A);nb=length(b);%当方程组行与列的维数不相等时，停止计算，并输出出错信息ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息ifm~=nberror('ThecolumnsofAmustbeequalthelengthofb!')return;end%开始计算，先赋初值flag='OK';x=zeros(n,1)