东南大学数值分析编程作业.docx
- 文档编号:11525192
- 上传时间:2023-03-12
- 格式:DOCX
- 页数:42
- 大小:543.50KB
东南大学数值分析编程作业.docx
《东南大学数值分析编程作业.docx》由会员分享,可在线阅读,更多相关《东南大学数值分析编程作业.docx(42页珍藏版)》请在冰豆网上搜索。
数值分析编程作业
姓名:
卢一鸣
学号:
151106
指导老师:
吴宏伟
目录
P20.第17题:
舍入误差与有效数 3
P56.第20题:
Newton迭代法 5
P126.第39题:
列主元Gauss消去法 13
P127.第40题:
逐次超松弛迭代法 16
P195.第37题:
3次样条插值函数 21
P256.第23题:
重积分的计算 23
P307.第23题:
常微分方程初值问题数值解 25
P346.第10题:
抛物方程Crank-Nicolson格式 30
注:
源程序放在对应的文件夹里,如第一题放在“P20.第17题:
舍入误差与有效数”文件夹里;
P20.第17题:
舍入误差与有效数
设
(1)编制按从小到大的顺序,计算的通用程序;
(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),'真实值']);
forN=[100100001000000]
w=sprintf('%s%d','N=',N);
disp([blanks(15),w]);
computation(N);
end
goon();
function[a]=computation(N)
%本函数计算三种SN值:
N按从小到大,N按从大到小,真实值
%第二个for循环N按从大到小顺序
%第二个for循环N按从大到小顺序
SN=0;
SN=single(SN);
%N从小到大的顺序
fori=2:
N
SN=SN+1/(i*i-1);
end
a=sprintf('%5.20f',SN);
disp([blanks(10),a]);
SN=0;
SN=single(SN);
%N从大到小的顺序
fori=N:
-1:
2
SN=SN+1/(i*i-1);
end
a=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]);
end
function[a]=goon()
%UNTITLED3Summaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
while1
disp('现在你可以尝试其他的N值,如想结束程序,请输入N<=1');
N=input('请输入N:
');
ifN<=1
return;
else
w=sprintf('%s%d','N=',N);
disp([blanks(15),w]);
computation(N);
end
end
(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从大到小变化时,可以避免大数吃小数,这在我们进行程序设计时,应尤其注意。
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;
newton('f','df',x0,delta);
%delta=search('f','df',delta)
function[y]=f(x)
%f(x)
y=x^3/3-x;
end
function[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);
while1
k=k+1;
x1=x0-feval('f',x0)/feval('df',x0);
err=abs(x1-x0);
x0=x1;
y=feval('f',x1);
if(err w=sprintf('%d\t%5.10f\t%e\t%5.10f',k,x0,err,y); disp(w); break end end function[delta]=search(f,df,delta) %寻找尽可能大的delta flag=1; k=1; x0=0; whileflag==1 delta=k*5e-4; x0=delta; k=k+1; m=0; flag1=1; whileflag1==1&&m<=10^3 x1=x0-feval('f',x0)/feval('df',x0); ifabs(x1-x0)<5e-4 flag1=0; end m=m+1; x0=x1; end ifflag1==1||abs(x0)>=5e-4 flag=0; end end end function[y]=df(x) %f'(x) y=x^2-1; end (2)①取容许误差为5e-4,通过search.m文件(执行命令delta=search('f','df',delta))找到尽可能大的。 ②。 区间上取-1000,-100,-50,-30,-10,-7,-5,-2,-1.5,-1.1,可见这些初始值都收敛于。 在区间即区间(-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.73205,(-δ,δ)区间收敛于0,(δ,1)区间类似于(-1,δ)区间,一部分收敛于1.73205,一部分收敛于-1.73205,(1,∞)收敛于1.73205。 通过本上机题,明白了对于多根方程,Newton法求方程根时,迭代序列收敛于某一个根有一定的区间限制,在一个区间上,可能会局部收敛于不同的根。 所以当我们求解这种方程时,最后先通过matlab的画图功能把曲线画出来,大概掌握零点的个数、位置以判断迭代出来的解是否正确和完整。 总的来说,初始值距离收敛值越近,迭代次数越小,但是当差距很远(如1000)时,迭代次数也不是太多,大概20次,所以newton迭代法只需较少的迭代次数就能得到收敛值。 P126.第39题: 列主元Gauss消去法 对于某电路的分析,归结为求解线性方程组,其中 (1)编制解n阶线性方程组的列主元Gauss消去法的通用程序; (2)用所编程序解线性方程组,并打印出解向量,保留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’表示计算成功 A,b [n,m]=size(A);nb=length(b); %当方程组行与列的维数不相等时,停止计算,并输出出错信息 ifn~=m error('TherowsandcolumnsofmatrixAmustbeequal! '); return; end%当方程组与右端项的维数不匹配时,停止计算,并输出出错信息 ifm~=nb error('ThecolumnsofAmustbeequalthelengthofb! ') return; end %开始计算,先赋初值flag='OK';x=zeros(n,1); fork=1: n-1 %选主元 max1=0; fori=k: n ifabs(A(i,k))>max1 max1=abs(A(i,k));r=i; end end ifmax1<1e-10 flag='failure';return; end %交换两行 ifr>k forj=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 %肖元过程 fori=k+1: n l=A(i,k)/A(k,k); forj=k+1: n A(i,j)=A(i,j)-l*A(k,j); end b(i)=b(i)-l*b(k); end end %回代过程 ifabs(A(n,n))<1e-10 flag='failure';return; end fork=n: -1: 1 forj=k+1: n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); end disp('x='); fprintf('%.5g\n',x) end function[A,b]=Aandb() %系数矩阵和右端向量 %Detailedexplanationgoeshere 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=[-15;27;-23;0;-20;12;-7;7;10]; %A=[2-46; %4-92; %1-13;]; %b=[3;5;4]; end (2)首先将R和V的内容输入到Aandb.m文件中。 输入命令main即得结果。 (3)本题编程过程中,由于是一段通用程序,要考虑到每一种可能,比如有计算成功的,有计算不成功的。 矩阵数据较多,在输出界面上最好能够再一次显示矩阵内容,方便用户查看矩阵输入是否正确。 要善于利用循环,掌握好每一个变量的意义,顺序,一个小地方出错则很可能导致整个程序运行不起来。 所以说编程是一项细活。 P127.第40题: 逐次超松弛迭代法 (1)编制解n阶线性方程组的SOR方法的通用程序(要求); (2)对于第39题中所给的线性方程组,取松弛因子,容许误差,打印松弛因子、迭代次数、最佳松弛因子及解向量。 答: (1)通用程序在SOR.m文件里。 程序执行前在Aandb.m文件中输入系数矩阵A和右端向量,执行main命令,即能获得结果。 %此函数为main函数 clearall [A,b]=Aandb(); delta=0.5e-5; max=10000;%如果迭代次数超过10000,失败 disp('wµü´ú´ÎÊý'); %[x,k,flag]=SOR(A,b,delta,1.96,max); fori=1: 99 w=i/50; [x,k,flag]=SOR(A,b,delta,w,max); a=sprintf('%1.2f\t%1.0f',w,k); disp(a); q(i)=w; c(i)=k; end [m,n]=min(c); [x,k,flag]=SOR(A,b,delta,q(n),max); a=sprintf('%s%1.2f\n%s','最佳松弛因子: ',q(n),'解向量'); disp(a); fprintf('%5.10g\n',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! ’表示迭代失败 ifnargin<5 max=10000; end ifnargin<4 w=1; end ifnargin<3 delta=1e-5; end n=length(A);k=0; x=zeros(n,1);y=zeros(n,1);flag='OK! '; while1 y=x; fori=1: n z=b(i); forj=1: n ifj~=i z=z-A(i,j)*x(j); end end ifabs(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 ifnorm(y-x,inf) break; end k=k+1; end function[A,b]=Aandb() %系数矩阵和右端向量 %Detailedexplanationgoeshere 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=[-15;27;-23;0;-20;12;-7;7;10]; %A=[430; %34-1; %0-14]; %b=[24;30;-24]; end (2)计算结果如下。 结果中迭代次数出现10000表示迭代了10000次还没能达到预定精度,故其不再迭代下去。 体会: 本程序将松弛因子取较小的步长,逐个求解,根据求解结果确定出最佳松弛因子,这种循环浪费计算机时间,但是貌似这也是唯一的方法。 也正是由于计算机,我们在求最佳松弛因子时,才会很方便。 ******************************************************************************* ************************************************************************************************************************************************************** 以上部分在第一次交作业时已交,以下部分是第二次作业******************************************************************************* ******************************************************************************* ************************************************************************************************************************************************************** ******************************************************************************* ************************************************************************************************************************************************************** ********************************************************************************************************************************************************************************************************************************************* ******************************************************************************* ************************************************************************************************************************************************************** ******************************************************************************* ************************************************************************************************************************************************************** ******************************************************************************* ************************************************************************************************************************************************************** ****************************
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 东南大学 数值 分析 编程 作业