数值分析课程课程设计

1、课 程 设 计 设计题目 数值分析 学生姓名 李飞吾 学 号 xxxxxxxx 专业班级 信息计xxxxx班 指导教师 设 计题 目共15题如下成绩课 程 设 计 主 要 内 容设计目的：通过不同题目的理解，进行算法分析。通过MATLAB软件进行编程对题目进行解决。个别题目设计验证，加深对数值分析的理解。函数的图像绘制的运用设计题目：题11 利用逆向递推的方法求解问题，通过条件终止地推题12 从某个初始值开始，利用递推公式进行积分估值题13 绘制Koch分形曲线，节点之间的关系与坐标变换题21 用高斯消元法的消元过程作矩阵分解，LU分解题22 矩阵分解方法求上题中A的逆矩阵，针对不同的b,而重

2、复利用已知的LU题23 验证希尔伯特矩阵的病态性,矩阵基本运算题31 用泰勒级数的有限项逼近正弦函数，由图像观察逼近效果题32 绘制飞机的降落曲线，线性方程组求解，与绘图题41 线性拟合的函数表达式的推导，使用了两种代码方法题51 用几种不同的方法求积分，观察数值积分的逼近效果题55 求空间曲线L弧长。求导后使用符号函数积分计算题61 用欧拉公式和四阶龙格-库塔法分别求解下列初值问题，代码搜索内容。题64 常微分方程的解，dsolve()函数使用题82 差分法解常微分方程边值问题，ode函数无能为力，Matlab中提供bvp解算器。 solinit=bvpinit(x,yinit,params

3、)sol= bvpsolver(odefun,bcfun,solinit,options) 题83 求解圆的半径，圆心。线性方程组解参数设计总结：（1） 算法是题目的解题核心，好的算法可以使计算更加精确 （题5.1）（2） 图形绘制在今后的课程设计，或者是论文中可以用到。（3） 无法解决的问题，需要请教室友，或者上网查阅。（4） MATLAB是一个很强大的软件，提供了很多内置的数学函数，直接进行解题。查阅资料时了解到一些MATLAB论坛。通过帖子阅读，了解到了MATLAB在科学计算方面，模拟，图形，视频等起到的作用。增加了对“计算科学“的理解。指 导 老 师 评 语建议：从学生的工作态度、工作

4、量、设计(论文的)创造性、学术性、使用性及书面表达能力等方面给出评价。 签名： 20 年 月 日 数值分析课程设计11 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？（15621）试分析椰子数目的变化规律，利用逆向递推的方法

5、求解这一问题解：算法分析：解该问题主要使用递推算法，关于椰子数目的变化规律可以设起初的椰子数为，第一至五次猴子在夜里藏椰子后，椰子的数目分别为再设最后每个人分得x个椰子,由题: （k=0,1,2,3,4）所以,利用逆向递推方法求解 （k=0,1,2,3,4）输出结果 ： 1023 15621MATLAB代码：n=input('n= ');n= 15621for x=1:n p=5*x+1; for k=1:5p=5*p/4+1; end结果分析：此题的解题思想是在迭代法中，判断p为整数时，输出与p if p=fix(p),breakendenddisp(x,p)12 设，（1）

6、从尽可能精确的近似值出发，利用递推公式：计算机从到的近似值；（2）从较粗糙的估计值出发，用递推公式：计算从到的近似值； 解：首先第一步，估计和的值：syms x n; int (x0/(5+x),0,1)ans=log(2)+log(3)-log(5)eval(ans)ans=0.1823则取为0.18syms x n;int(x30/(5+x),0,1)ans =931322574615478515625*log(2)+931322574615478515625*log(3)-931322574615478515625*log(5)-79095966183067699902965545527

7、073/465817912560eval(ans)ans = 0MATLAB中小数点后保留四位，由上面计算知道积分的值不为了零。 所以的取值为0.00001-0.0001MATLAB代码：i=input('i=');i=0.18;>> if i>=0.1&&i<=0.2 for n=1:1:20 i=-5*i+1/n endelseif i>0&&i<=0.0001 for n=30:-1:2 i=(-1/5)*i+1/(5*n) endend i = 1.1336e+005i = -5.6679e+005i

8、= 2.8339e+006i = -1.4170e+007i = 7.0848e+007i = -3.5424e+008i = 1.7712e+009i = -8.8560e+009i = 4.4280e+010i = -2.2140e+011同理输入积分初始值i=0时可以得 i=0.0884结果分析：第二种方法所得的结果相对来说比较精确一些，也比较可靠因为第一种方法每一迭代都将最初的误差放大了五倍，使得最终的误差越来越大；而第一种方法经过每一次迭代都将误差缩小为初始误差的五分之一，使得最终的误差越来越小，因此相对来说比较可靠，性能较好。13 绘制Koch分形曲线问题描述：从一条直线段开始，将

9、线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形（图1-4）；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形（图1-5），这时，图形中共有17个结点。这种迭代继续进行下去可以形成Koch分形曲线。问题分析：考虑由直线段（2个点）产生第一个图形（5个点）的过程，设和分别为原始直线段的两个端点。现在需要在直线段的中间依次插入三个点产生第一次迭代的图形（图1-4）。显然，位于点右端直线段的三分之一处，点绕旋转60度（逆时针方向）而得到的，故可以处理为向量经正交变换而得到向量，形成算法如下：（1）；（2）；（3）；在

10、算法的第三步中，A为正交矩阵。这一算法将根据初始数据（和点的坐标），产生图1-4中5个结点的坐标。这5个结点的坐标数组，组成一个5×2矩阵。这一矩阵的第一行为为的坐标，第二行为的坐标，第二行为的坐标第五行为的坐标。矩阵的第一列元素分别为5个结点的x坐标 ，第二列元素分别为5个结点的y坐标。问题思考与实验：（1）考虑在Koch分形曲线的形成过程中结点数目的变化规律。设第k次迭代产生结点数为，第迭代产生结点数为，试写出和之间的递推关系式；（2）参考问题分析中的算法，考虑图1-4到图1-5的过程，即由第一次迭代的5个结点的结点坐标数组，产生第二次迭代的17个结点的结点坐标数组的算法；（3）

11、考虑由第k次迭代的个结点的结点坐标数组，产生第次迭代的个结点的结点坐标数组的算法；（4）设计算法用计算机绘制出如下的Koch分形曲线（图1-6）解：(1) (2)(3)算法及(4)代码分析：p=0 0;10 0; %P为初始两个点的坐标,第一列为x坐标,第二列为y坐标n=2; %n为结点数A=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); %旋转矩阵for k=1:4 d=diff(p)/3; %diff计算相邻两个点的坐标之差,得到相邻两点确定的向量 %则d就计算出每个向量长度的三分之一,与题中将线段三等分对应 m=4*n-3; %迭代公式 q=p(1:n

12、-1,:); %以原点为起点,前n-1个点的坐标为终点形成向量 p(5:4:m,:)=p(2:n,:); %迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标 p(2:4:m,:)=q+d; %用向量方法计算迭代后处于4k+2位置上的点的坐标 p(3:4:m,:)=q+d+d*A' %用向量方法计算迭代后处于4k+3位置上的点的坐标 p(4:4:m,:)=q+2*d; %用向量方法计算迭代后处于4k位置上的点的坐标 n=m; %迭代后新的结点数目end plot(p(:,1),p(:,2) %绘出每相邻两个点的连线axis(0 10 0 3.5)21 用高斯消元法的消元过程作矩阵分解

13、。设消元过程可将矩阵A化为上三角矩阵U，试求出消元过程所用的乘数、并以如下格式构造下三角矩阵L和上三角矩阵U验证：矩阵A可以分解为L和U的乘积，即A=LU。矩阵LU分解MATLAB代码:function hl=zhjLU(A)n,n=size(A);RA=rank(A);if RA=n disp('因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下：'); RA,hl=det(A); returnendif RA=n for p=1:n h(p)=det(A(1:p,1:p); end hl=h(1:n); for i=1:n if h(1,i)=0 disp

14、('因为A的各阶主子式等不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'); RA,hl return end end if h(1,i)=0 disp('因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl如下：'); for j=1:n U(1,j)=A(1,j); end for k=2:n for i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if i>j L(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1); L(i

15、,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k)/U(k,k); else U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end RA,hl,U,L· endend以上上代码保存为M文件，并在命令窗口输入A=20 2 3;1 8 1; 2 -3 15;b=0 0 0'h=zhjLU(A)程序运行结果：L = U= 1.0000 0 0 20.0000 2.0000 3.0000 0.0500 1.0000 0 0 7.9000 0.8500 0.1000 -0.4051 1.0000 0 0 15.0443

16、实验验证： 可以直接使用MATLAB内置LU分解A=20 2 3;1 8 1; 2 -3 15;L U=lu(A) 输出结果与上程序输出结果一致。22 用矩阵分解方法求上题中A的逆矩阵。记 分别求解方程组 由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。对第一个方程组，由于A=LU，所以先求解下三角方程组,再求解上三角方程组，则可得逆矩阵的第一列列向量；类似可解第二、第三方程组，得逆矩阵的第二列列向量的第三列列向量。由三个列向量拼装可得逆矩阵。解:MATLAB代码如下：b1=1;0;0; b2=0;1;0; b3=0;1;1;A=20,2,3;1,8,1;2,-3,15;L=1,0,0

17、;0.05,1,0;0.1,-0.4051,1;U=20 2 3;0 7.9 0.85;0 0 15.0443;Y1=Lb1 X1=UY1 Y2=Lb2X2=UY2Y3=Lb3X3=UY3Y1 = 1.0000 -0.0500 -0.1203X1 = 0.0517 -0.0055 -0.0080Y2 =实验验证： X1 X2 X3ans = 0.0517 -0.0164 -0.0257 -0.0055 0.1237 0.1165 -0.0080 0.0269 0.0934而：inv（A）=X1 X2 X3 得证 0 1.0000 0.4051X2 = -0.0164 0.1237 0.0269

18、Y3 = 0 1.0000 1.4051X3 = -0.0257 0.11650.0934- 19 -23 验证希尔伯特矩阵的病态性：对于三阶矩阵取右端向量，验证：（1）向量是方程组的准确解；（2）取右端向量b的三位有效数字得，求方程组的准确解，并与X的数据作比较 。说明矩阵的病态性。解：(1)H=1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5;X=1;1;1;b=H*Xb = 1.8333 1.08330.7833 与题中相同(2) 先求出解X，与数据作比较 。H=1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5;b=1.83;1.08;0.783;X=

19、HbX = 1.0800 0.54001.4400与相差较大，矩阵为病态矩阵31 用泰勒级数的有限项逼近正弦函数用计算机绘出上面四个函数的图形。解：MATLAB代码如下(1) syms x;taylor(sin(x)x=0:0.01*pi:piplot(x,sin(x) (2) syms x;taylor(x)x=0:0.01*pi:pi/2plot(x) (3) syms x;taylor(x-x3/6)fplot('x-x3/6',0 pi/2) (4) syms x;taylor(x-x3/6+x5/120)fplot('x-x3/6+x5/120',0

20、pi/2) 结果图形右： x=0:0.1:pi;y=sin(x);plot(x,y,'-sk');hold onx=0:0.1:pi/2;y=x;plot(x,y,'-b*')hold onfplot('x-x.3/6',0,pi/2,0,2,2e-3,'-gx')hold onfplot('x-x.3/6+x.5/120',0,pi/2,0,2,2e-3,'-r.')hold offlegend('sin(x)','x','x-x3/6','

21、x-3/6+x5/120',2)xlabel('x')ylabel('y')title('Taylor approximation')结果分析：图中红色点线为正弦曲线，蓝色的星线为一阶泰勒逼近，绿色叉线为二阶泰勒逼近，黑色正方形线为三阶泰勒逼近。可见三阶泰勒逼近效果最好，泰勒级数越高，逼近效果越好。32 绘制飞机的降落曲线一架飞机飞临北京国际机场上空时，其水平速度为540km/h，飞行高度为1 000m。飞机从距机场指挥塔的横向距离12 000m处开始降落。根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。建立直角坐标系，设飞机着陆点

22、为原点O，降落的飞机为动点，则表示飞机距指挥塔的距离，表示飞机的飞行高度，降落曲线为该函数满足条件：（1）试利用满足的条件确定三次多项式中的四个系数；（2）用所求出的三次多项式函数绘制出飞机降落曲线。function s=f(x1,y1,x2,y2,x3,y3,x4,y4)format longa1=1 ,x1,x12,x13;a2=1 ,x2,x22,x23;a3=0,1,2*x3,3*x32;a4=0,1,2*x4,3*x42;a=a1; a2;a3;a4;b=y1; y2;y3;y4;s=ab;x=-12000:250:0;y=s(3)*x.2-s(4)*x.3plot(x,y,'

23、;-d')xlabel('x')ylabel('y')以上为M文件内容，在命令窗口键入f(0,0,12000,1000,0,0,12000,0)运行结果： ans为多项式系数ans = 1.0e-004 * 0 0 0.20833333333333 -0.0000115740740741 曾任英特尔公司董事长的摩尔先生早在1965年时，就观察到一件很有趣的现象：集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。因而发表论文，提出了大大有名的摩尔定律（Moores Law），并预测未来这种增长仍会延续下去。下面数据中，第二行数据为晶片

24、上晶体数目在不同年代与1959年时数目比较的倍数。这些数据是推出摩尔定律的依据：年代19591962196319641965增加倍数13456试从表中数据出发，推导线性拟合的函数表达式。解：解法一MATLAB代码：x=1959,1962,1963,1964,1965;y=1,3,4,5,6;p1=polyfit(x,y,1)y1=polyval(p1,x)plot(x,y1,'-',x,y,'r*')xlabel('x'),ylabel('y');运行结果：p1 = 1.0e+003 * 0.0008 -1.6255y1 =0.8

25、113 3.3019 4.1321 4.9623 5.7925线性拟合的函数表达式： Y=0.8302x-1.6255e+003解法二：运行结果：xs = 1.0e+003 * -1.625528301891238 0.000830188679248x=1959 1962 1963 1964 1965;y=1;3;4;5;6;for i=1:length(x) for j=1:2 A(i,j)=x(i)(j-1); endendL,U=lu(A'*A);xs=U(L(A'*y)从而年代y与增加倍数x之间的关系为：y=-1625.528301891238+0.8301886792

26、48x51 用几种不同的方法求积分的值。（1）牛顿-莱布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）复合梯形公式。解：syms xi1=int(4/(1+x2),x,0,1)运行结果：i1 =pii2 = 3i3 = 3.1333i4 = 3.1416a=0;b=1;h=b-a;i2=(4/(1+a2)+4/(1+b2)/2i3=h/6*(4/(1+a2)+4*4/(1+(a+b/2)2)+4/(1+b2)M=100;h=(b-a)/M;i4=0;for k=1:(M-1) x=a+h*k; i4=i4+4/(1+x2);endi4=h*(4/(1+a2)+4/(1+b2)/2+h*i4

27、结果分析：牛顿莱布尼兹公式得到精确结果，采用梯形公式得到的结果比采用Simpson公式的精确度要低，采用复化梯形公式在步长取得越来越小的状态下可以提高精度。55 求空间曲线L：弧长公式为解：运行程序为：syms t;x=diff(cos(t);y=diff(sin(t);z=diff(2-cos(t)-sin(t);y=int(x2+y2+z2)0.5,t,0,2*pi); %符号函数积分digits(14);i=vpa(y)运行结果：i = 8.737752570989461 用欧拉公式和四阶龙格-库塔法分别求解下列初值问题；解：（1）<1>欧拉公式：function t,x=E

28、uler(fun,t0,tt,x0,N)h=(tt-t0)/N;t=t0+0:N'*h;x(1,:)=x0'运行结果：ans = 0 1.0000 0.0500 1.0450 0.1000 1.0878 0.1500 1.1285 0.2000 1.1676 0.2500 1.2051 0.3000 1.2413 0.3500 1.2762 0.4000 1.3100 0.4500 1.3427 0.5000 1.3745 0.5500 1.4055 0.6000 1.4356 0.6500 1.4649 0.7000 1.4936 0.7500 1.5216 0.8000 1

29、.5490 0.8500 1.5758 0.9000 1.6021 0.9500 1.6278 1.0000 1.6531for k=1:N f=feval(fun,t(k),x(k,:); f=f' x(k+1,:)=x(k,:)+h*f;end以'Euler.m'保存function f=Euler_fun(t,x)f=0.9*x./(1+2*t);以'Euler_fun.m'保存function main_Eulert,x=Euler('Euler_fun',0,1,1,20);fh=dsolve('Dx=0.9*x/(1+

30、2*t)','x(0)=1');for k=1:8 ft(k)=t(k); fx(k)=subs(fh,ft(k);endt,x以'main_Euler.m'保存输入：main_Euler运行结果：ans = 0 1.0000 0.1250 1.0954 0.2500 1.1807 0.3750 1.2585 0.5000 1.3306 0.6250 1.3981 0.7500 1.4617 0.8750 1.5223 1.0000 1.5801<2>四阶龙格-库塔法function R=rk4(f,a,b,ya,N)%y'=f(x,

31、y)%a,b为左右端点%N为迭代步长%h为步长%ya为初值h=(b-a)/N;T=zeros(1,N+1);Y=zeros(1,N+1);T=a:h:b;Y(1)=ya;for j=1:N k1=h*feval(f,T(j),Y(j); k2=h*feval(f,T(j)+h/2,Y(j)+k1/2); k3=h*feval(f,T(j)+h/2,Y(j)+k2/2); k4=h*feval(f,Y(j)+h,Y(j)+k3); Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endans=T' Y' 以'rk4.m'保存function z=f

32、(x,y)z=0.9*y./(1+2*x);以'f.m'保存输入：rk4('f',0,1,1,8)64 列出函数在区间0，e上的函数值表并作出它的图象。其中，是初值问题 的解。解v=dsolve('Dv*log(x)=2*x','v(0)=0','x');f=(1-log(v)*v f =-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x) ezplot('f')输出结果：f =-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x)7.4

33、一个10次项式的系数为1 a1 a2 a9 a10=1 55 1320 181 50 157 773 902 055 341 693 0 -840 950 0 127 535 76 -106 286 40 632 880 0试用多项式的求根指令roots求出该10次方程的10个根，然后修改9次项的系数-55为-56，得新的10次方程，求解新的方程，观察根的变化是否很显著。解： p=1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;>> roots(p)ans = 10.6051

34、+ 1.0127i 10.6051 - 1.0127i 8.5850 + 2.7898i 8.5850 - 2.7898i 5.5000 + 3.5058i 5.5000 - 3.5058i 2.4150 + 2.7898i 2.4150 - 2.7898i 0.3949 + 1.0127i 0.3949 - 1.0127ip=1 -56 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;>> roots(p)ans = 21.7335 7.3501 + 7.7973i 7.3501 - 7.7973i 4.3589 + 3.3285i 4.3589 - 3.3285i 5.1831 2.4378 + 2.7974i 2.4378 - 2