数值分析课程设计报告.doc

本科课程设计（报告）数值分析课程设计设计题一 设计实验验证Hilbert矩阵的病态性。 （提示：先取x（比如x=（1,1.1）计算出b=Hx,然后通过列主元GAUSS消去法求解Hx=b,得到近似解,比较之，并研究随着n增大，解的误差变化情况，得出结论 ）设计思路： （1）取定初值X=（1，1，1，），计算出b=Hx,然后通过列主元GAUSS消去法求解Hx=b,得到近似解B=（b1,b2,b3）； （2）再把近似解B，代入方程b=Hx,，解出X1，比较X与X1，判断H矩阵的病态性； （3）比较各个H矩阵的病态性，研究它们的误差变化情况。 算法步骤： （1）取一个三阶H矩阵，按照设计思路的方法，判断出三阶H矩阵的病态性； （2）取一个四、五阶H矩阵，同理，判断出四阶H矩阵的病态性。 程序清单：Xa=1;1;1;Ha=1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5;Ba=Ha*XaBa=1.83;1.08;0.78;Xa=HaBaXb=1;1;1;1;Hb=1,1/2,1/3,1/4;1/2,1/3,1/4,1/5;1/3,1/4,1/5,1/6;1/4,1/5,1/6,1/7;Bb=Hb*XbBb=2.08;1.28;0.95;0.76;Xb=HbBbXc=1;1;1;1;1;Hc=1,1/2,1/3,1/4,1/5;1/2,1/3,1/4,1/5,1/6;1/3,1/4,1/5,1/6,1/7;1/4,1/5,1/6,1/7,1/8;1/5,1/6,1/7,1/8,1/9;Bc=Hc*XcBc=2.28;1.45;1.09;0.89;0.75;Xc=HcBc 运行过程与输出结果截图： 结果分析：代入近似解Ba,Bb,Bc而得到的X 如下：Xa = 0.99000000000000 1.08000000000001 0.89999999999999Xb = 1.27999999999997 -1.79999999999961 7.19999999999903 -2.79999999999936Xc = 1.0e+002 * -0.06999999999975 1.48199999999542 -6.25799999998035 9.37999999997045 -4.53599999998560由此分析比较初始的Xa,Xb,Xc可以看出，Hilbert矩阵是一个病态矩阵，且随着n的增大，其误差变化加大。即，Hilbert矩阵是一个随n增大而误差偏离更大的病态矩阵。设计题二 1225年，达。芬奇研究了方程并得到它的一个近似根。没有人知道他用什么方法得到它。设计两种方法去计算，并比较这两种方法。设计思路与算法步骤：不动点迭代法：（1） 先把方程变换为x=20/(x.2+2*x+10)，并编译函数M文件（2） 利用不动点迭代法，确定迭代次数，进行运算，不妨设迭代次数为20，由题目可设迭代初始值X0=1牛顿迭代法：（1） 用牛顿迭代法，得到迭代方程 x1=x0f(x0)/f(x0)，编译函数M文件（2）取适当的，和迭代次数n,进行迭代，直到符合条件时，停止运算。（3）确定迭代初始值x0=1，设精确为0.000000001.程序清单：1、 函数文件：原函数;function y = fun(x)y=x.3+2*x.2+10*x-20;迭代函数：function y=fun1(x)y=20/(x.2+2*x+10);导函数：function y =dfun(x)y=3*x.2+4*x+10; 2、不动点迭代的M文件：function k,pcha,xpcha,xk=diedai(x0,k)% 输入的量-x0是初始值,k是迭代次数x(1)=x0; for i=1:k x(i+1)=fun1(x(i);%程序中调用的fun1.m为函数y=(x) pcha= abs(x(i+1)-x(i); xpcha=pcha/( abs(x(i+1)+eps);%偏差计算 i=i+1;xk=x(i);(i-1) pcha xpcha xk end if (pcha 1)&(xpcha0.5)&(k3) %迭代收敛性判断 disp(请用户注意：此迭代序列发散，请重新输入新的迭代公式) return; end if (pcha 0.001)&(xpcha3) disp(祝贺您！此迭代序列收敛，且收敛速度较快) return; endp=(i-1) pcha xpcha xk;牛顿迭代的M文件：function k,xk,yk,pcha,xpcha=newton(x0,tx,fx,n) %k为迭代初始值，xpcha，pcha是误差范围，迭代终止条件，n控制迭代次数x(1)=x0; for i=1: n x(i+1)=x(i)-fun(x(i)/(dfun(x(i)+eps); %牛顿迭代公式与偏差计算表达式pcha=abs(x(i+1)-x(i); xpcha= pcha/( abs(x(i+1)+eps); i=i+1;xk=x(i);yk=fun(x(i); (i-1) xk yk pcha xpcha if (abs(yk)fx)&(pchatx)|(xpchan k=i-1; xk=x(i);(i-1) xk yk pcha xpcha return;end (i-1),xk,yk,pcha,xpcha运行过程与输出结果截图： %在命令窗口中输入以下语句，可得以下结果:k,pcha,xpcha,xk=diedai(1,20) k,xk,yk,pcha,xpcha=newton(1,0.000000001, 0.000000001,1000)不动点迭代牛顿迭代结果分析：迭代法：运行结果k = 20pcha = 1.077823423845103e-007xpcha = 7.874174939314717e-008xk = 1.36880807468942容易看出，进行20次迭代之后，得到的结果x = 1.36880807468942，而且与题目所给的，十分接近。牛顿迭代法：运行结果k = 5xk = 1.36880810782137yk = 0pcha = 1.7863568394251e-015xpcha = 1.2987398576591e-015可知，利用牛顿迭代法进行五次迭代后，即可得到满足精度的结果，与题目给的答案接近。 总而言之，不动点迭代法与牛顿迭代法这两种方法，可见，牛顿迭代法，运算次数少，且误差小于不动点迭代的运算误差，故牛顿迭代法比较优越。设计题四 地球卫星飞行的轨道是一个椭圆，椭圆的周长计算公式为 其中a为椭圆长半轴，c是椭圆中心与地球中心的距离。令H,h分别为远地点，和近地点距离。R=6371（km）为地球半径。则a=(2R+H+h)/2, c=(H-h)/2.我国第一颗人造卫星H=3484km,h=439km.（a=(2*6371+3484+439)/2=8332.5, c=(3484-439）/2)=1522.5试通过复化梯形公式，复化辛普森公式计算卫星轨道的周长，且研究并比较这两种的方法的收敛性。设计思路与算法步骤：(1)、编译被积函数M文件，被积函数表达式为：y=8332.5*sqrt(1-(1522.5/8332.5)2*sin(x);（2）、根据复化梯形公式， 编译算法。不妨设递归次数为15次。 （由于复化梯形公式，的两个端点值，所使用的公式不同，故要分开编写循环体）（3）、根据复化辛普森公式，编译算法，不妨设递归次数为15次。程序清单：被积函数：function y=fun4(x)y=8332.5*sqrt(1-(1522.5/8332.5)2*sin(x);复化梯形公式;function T=rctrap(fun4,a,b,m) %fun4表示被积函数，a，b是被积函数上下限，m是递归次数n=1;h=b-a; T=zeros(1,m+1); x=a; %赋值计算T(1)=h*(feval(fun4,a)+feval(fun4,b)/2; %初次调用梯形公式for i=1:m %复化梯形公式主体 h=h/2; n=2*n; s=0; for k=1:n/2 x=a+h*(2*k-1); s=s+feval(fun4,x);endT(i+1)=T(i)/2+h*s; %最后一次调用复化梯形公式endT=T(1:m);复化辛普森公式：function y,n,yiter=SimpsonInt(fun,a,b,epsilon)%Function : Simpson积分%Argument : fun是被积函数；a,b是积分上下限；epsilon是误差容限%Return : y为积分结果；将a,b平均n等分；yiter迭代的积分值%if (nargin4)epsilon=1e-6;endif (nargin3)warning(usage:,mfilename,(fun,a,b,epsilon);endn=2;h=(b-a)/4;y1=feval(fun,a)+feval(fun,b);y2=feval(fun,(a+b)/2);y0=(b-a)/6*(y1+4*y2);yiter=y0;while 1y3=sum(feval(fun,a+(1:2:2*n-1)*h);y=h/3*(y1+2*y2+4*y3);if abs(y-y0)15*epsilonbreak; endn=n+n;h=h/2;y2=y2+y3;y0=y;yiter=yiter,y0;end运行过程与输出结果截图： 在命令窗口中输入：Q1=rctrap(fun4,0,pi/2,15)vpa(Q1,7)y,n,yiter=SimpsonInt(fun4,0,pi/2,0.01)vpa(y,7)可得以下截图结果分析：复化梯形公式vpa(Q1,7)ans = 12978.49, 12955.87, 12950.43, 12949.09, 12948.75, 12948.67, 12948.65, 12948.64, 12948.64, 12948.64, 12948.64, 12948.64, 12948.64, 12948.64, 12948.64复化辛普森公式vpa(y,7)ans =12948.64可见，利用复化梯形公式计算得到的周长为12948.64, 而利用复化辛普森公式计算得到的周长为12948.64时，仅递归了4次。由此可见，复化辛普森公式的收敛速度比复化梯形公式快。设计题五给定单摆方程初值问题其中g=9.8,l=25.取初始偏离角度取初始偏离角度其精确解为。分别对上述两种情况按照下列方法求出其数值解，比较各方法的优缺点，并将计算结果与精确解做比较（列表、画图）。（方案I）欧拉法，步长h = 0.025, h = 0.1；（方案II）改进的欧拉法，步长h = 0.05, h = 0.1；（方案III）四阶经典龙格库塔法，步长h = 0.1。设计思路与算法步骤： 对于高阶常微分方程的解法：1、把二阶微分方程化成标准型dt/dx=-g/l*sinx t(0)=0 x(0)= x=; t =x2、把函数编译成矩阵表达式3、解出精确解4、利用欧拉法，改进的欧拉法，四阶龙格库塔法，求解5、在相同的步长内，比较三种方法和精确解的误差大小 程序清单：函数M文件:function f=dy(x,y)f=y(2);-9.8*sin(y(1)/25;(1)精确解：y=0.1745*cos(9.8*x/25)plot(x,y,-);xlabel(x), ylabel(y)，（2） 欧拉法：function x,y=Euler(fun,x0,xt,y0,n,h);%xt是右端点值x(1)=x0;y=y0; %赋初值for i=1:n x(i+1)=x(i)+h; y(:,i+1)=y(:,i)+h*feval(fun,x(i),y(:,i); %欧拉公式endplot(x,y,-) %作图xlabel(x), ylabel(y), （3）改进的欧拉公式：function x,y=Impeuler(fun,x0,xt,y0,n,h);%定义改进的欧拉函数x(1)=x0;y=y0; %赋值for i=1:nx(i+1)=x(i)+h; %改进欧拉公式y1=y(:,i)+h*feval(fun,x(i),y(:,i);y2=y(:,i)+h*feval(fun,x(i+1),y1);y(:,i+1)=(y1+y2)/2;endplot(x,y,-) %作图xlabel(x), ylabel(y) (4) 四阶经典龙格库塔法：function x,y=RK4(fun,y0,h,a,n); x(1)=a; %赋值 y(:,1)=y0; for i=1:n %四阶经典龙格库塔法公式 x(i+1)=x(i)+h; k1=feval(fun,x(i),y(:,i); k2=feval(fun,x(i)+h/2,y(:,i)+h*k1/2); k3=feval(fun,x(i)+h/2,y(:,i)+h*k2/2); k4=feval(fun,x(i)+h,y(:,i)+h*k2); y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;endplot(x,y,-) %作图xlabel(x), ylabel(y), 运行过程与输出结果截图： %在命令窗口中输入，可得以下结果 x,y=euler(dy,0,1,0.1745;0,10,0.1) x,y=euler(dy,0,1,0.1745;0,40,0.025) x,y=euler(dy,0,1,0.5236;0,40,0.025) x,y=euler(dy,0,1,0.5236;0,10,0.1) x,y=Impeuler(dy,0,1,0.1745;0,20,0.05) x,y= Impeuler(dy,0,1,0.1745;0,10,0.1) x,y= Impeuler(dy,0,1,0.5236;0,20,0.05) x,y= Impeuler(dy,0,1,0.5236;0,10,0.1) x,y=RK4(dy,0.5236;0,0.1,0,10) x,y=RK4(dy,0.1745;0,0.1,0,10) y=0.1745*cos(sqrt(9.8/25)*x) y=0.5236*cos(sqrt(9.8/25)*x) 截图：1、 精确解：（方案I）欧拉法：（1） 步长h = 0.025， （2）h=0.1 (3) h=0.025 (4) h=0.1 （方案II）改进的欧拉法：（1） 步长h = 0.05, （2） h=0.1, (1) h=0.1, (2) h=0.05, （方案III）四阶经典龙格库塔法:(1) (2) 结果分析：1、当，h=0.1时，从上诉的程序运行过程和结果可以看出，利用改进的欧拉法和经典四阶龙格库塔公式解出的值，比较接近精确解，而欧拉法比较远离精确解。且经典四阶龙格库塔公式收敛相对快一点。2、当，h=0.1时， 得出的结论，与时，相同。总而言之，利用四阶龙格库塔方法所得的结果比较具有优越性，且收敛快。而欧拉算法，比较简单，运算比较快，改进的欧拉法处于两者之间心得体会这次课程设计的时间不算长，可是我从中学习到了很多东西。由于有选择数学软件这个课程，所以对于matlab这个软件的认识还是有一定的基础。不过，由于实际操作比较小，所以处理起这个课程设计还是比较吃力。从一开始，我对于matlab的一无所知，到现在可以利用这个软件对一些常用的算法，加于研究和设计，途中经历了很多困难。比如，如何使用for语句，矩阵的运算操作等等，都不太清楚。具体的，例如设计题二。虽然利用数值分析上所学到的知识，可以从理论上解决这个问题。但是要利用matlab来处理的话，我顿时觉得无从下手，因为，要设计一个算法，毕竟涉