数值分析课程设计.docx
- 文档编号:8613247
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:25
- 大小:188.88KB
数值分析课程设计.docx
《数值分析课程设计.docx》由会员分享,可在线阅读,更多相关《数值分析课程设计.docx(25页珍藏版)》请在冰豆网上搜索。
数值分析课程设计
数值分析
课程设计报告
专业:
班级:
学号:
姓名:
1.1水手、猴子和椰子问题:
五个水手带了一只猴子来到南太平洋的一个荒岛上,发现那里有一大堆椰子。
由于旅途的颠簸,大家都很疲惫,很快就入睡了。
第一个水手醒来后,把椰子平分成五堆,将多余的一只给了猴子,他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来,和第一个水手一样,把椰子分成五堆,恰多一只猴子,私藏一堆,再去入睡,天亮以后,大家把余下的椰子重新等分成五堆,每人分一堆,正好余一只再给猴子,试问原先共有几只椰子?
试分析椰子数目的变化规律,利用逆向递推的方法求解这一问题(15621)。
分析:
定义数组x()用于存储每天的椰子数
通过给最后一天的椰子数x(7)附初值根据公式“x(n)=5*x(n+1)/4+1”逆向递推前一天的椰子数
如果推出某一天的椰子数不为整数(即:
mod(5*x(n+1),4)~=0),说明初值x(7)不正确,则进行x(7)=x(7)+1运算来改变初值,直到符合要求为止。
最终一次性输x(),6次分配的前后值。
程序代码:
functionyezi1(m)
x(7)=0;
whilex(7) x(7)=m; x(6)=5*x(7)+1; n=5; while(mod(5*x(n+1),4)==0)&(n>0) x(n)=5*x(n+1)/4+1; n=n-1; end ifn~=0 m=x(7)+1; end end disp(x) 指令窗口输入: >>m=1;yezi1(m) 回车输出: 156211249699967996639651161023 本程序,首先赋值m=1(m为最后分完时每堆椰子数量),程序自身会测试m的可行性,如不满足要求,程序自动增大m值,直到出现第一个结果,此时m=1023,总椰子数为15261. 如果输入m初值>1023,程序会自动找到下一个符合值2047,………… 2.1用高斯消元法的消元过程作矩阵分解。 设 消元过程可将矩阵A化为上三角矩阵U,试求出消元过程所用的乘数 、 、 并以如下格式构造下三角矩阵L和上三角矩阵U 验证: 矩阵A可以分解为L和U的乘积,即A=LU。 分析: 首先用矩阵A的各阶主子式值判断矩阵能否进行LU分解,在能够分解的前提下,使用高斯消去法,将将消去时的两行间的倍数m附入L(m=A(k,p)/A(p,p);L(p,p)=1;L(k,p)=m;),再将消去后的A附入U,即得L,U。 程序代码: function[L,U]gaussLU(A) [nn]=size(A);RA=rank(A); ifRA~=n disp('因为A的n阶行列式hl等于零,所以A不能进行LU分解。 '); hl=det(A); return end ifRA==n forp=1: n h(p)=det(A(1: p,1: p)); end hl=h(1: n); fori=1: n ifh(1,i)==0 disp('因为A的r阶主子式等于零,所以A不能进行LU分解。 ') return end end ifh(1,i)~=0 disp('因为A的各阶主子式都不等于零,所以A能进行LU分解。 L与U的值如下: ') L=zeros(n,n);L(n,n)=1; U=zeros(n,n); forp=1: n-1 fork=p+1: n m=A(k,p)/A(p,p); L(p,p)=1;L(k,p)=m; A(k,p: n)=A(k,p: n)-m*A(p,p: n); end end U=A(1: n,1: n); end end 指令窗口输入: A=[2023;181;235];[L,U]=gaussLU(A) 程序输出: 因为A的各阶主子式都不等于零,所以A能进行LU分解。 L与U的值如下: L= 1.000000 0.05001.00000 0.10000.35441.0000 U= 20.00002.00003.0000 07.90000.8500 004.3987 2.2用矩阵分解方法求上题中A的逆矩阵。 记 分别求解方程组 由于三个方程组系数矩阵相同,可以将分解后的矩阵重复使用。 对第一个方程组,由于A=LU,所以先求解下三角方程组 再求解上三角方程组 ,则可得逆矩阵的第一列列向量;类似可解第二、第三方程组,得逆矩阵的第二列列向量的第三列列向量。 由三个列向量拼装可得逆矩阵 。 分析: 我觉得题目中已经分析的很全面了,此题没有必要再多此一举了。 程序代码: A=[20,2,3;1,8,1;2,-3,15]; [L,U]=gaussLU(A); b1=[1;0;0];b2=[0;1;0];b3=[0;0;1]; Y1=L\b1; Y2=L\b2; Y3=L\b3; X1=U\Y1; X2=U\Y2; X3=U\Y3; X=[X1,X2,X3]%X 输出结果: X= 0.0517-0.0164-0.0093 -0.00550.1237-0.0072 -0.00800.02690.0665 实验三 3.1用泰勒级数的有限项逼近正弦函数 用计算机绘出上面四个函数的图形。 分析: 直接设定变量范围、函数,然后作图即可,注意图像可能会重叠,最好用不同线条,区别大一点的颜色绘制。 程序代码: x=0: pi/100: pi; x1=0: pi/100: pi/2; y0=sin(x); y1=x1; y2=x1-x1.^3/6; y3=x1-x1.^3/6+x1.^5/120; plot(x,y0,': k',x1,y1,'-g',x1,y2,'-b',x1,y3,'-r') 输出图像: 3.2绘制飞机的降落曲线 一架飞机飞临北京国际机场上空时,其水平速度为540km/h,飞行高度为1000m。 飞机从距机场指挥塔的横向距离12000m处开始降落。 根据经验,一架水平飞行的飞机其降落曲线是一条三次曲线。 建立直角坐标系,设飞机着陆点为原点O,降落的飞机为动点 ,则 表示飞机距指挥塔的距离, 表示飞机的飞行高度,降落曲线为 该函数满足条件: (1)试利用 满足的条件确定三次多项式中的四个系数; (2)用所求出的三次多项式函数绘制出飞机降落曲线。 分析: 通过y(x)满足的条件带入降落曲线方程,即可得到一个四元一次方程组。 将方程组视为矩阵。 由于前面的题目用过了高斯消去法,本体继续沿用此方法,通过高斯消去法,求解方程组,得到a0,a1,a2,a3的值,带回y中,运用plot函数,输出x,y的关系曲线。 程序代码: functiona=jiangl(x1,x2,x3,x4,b) A=[1x1x1^2x1^3;1x2x2^2x2^3;012*x33*x3*x3;012*x43*x4*x4]; B=[Ab]; a=zeros(4,1);C=zeros(1,5); forp=1: 3 [Y,j]=max(abs(B(p: 4,p)));C=B(p,: ); B(p,: )=B(j+p-1,: );B(j+p-1,: )=C; fork=p+1: 4 m=B(k,p)/B(p,p); B(k,p: 5)=B(k,p: 5)-m*B(p,p: 5); end end b=B(1: 4,5);A=B(1: 4,1: 4);a(4)=b(4)/A(4,4); forq=3: -1: 1 a(q)=(b(q)-sum(A(q,q+1: 4)*a(q+1: 4)))/A(q,q); end 在指令窗口输入: >>x1=0;x2=12000;x3=0;x4=12000;b=[0;1000;0;0];a=jiangl(x1,x2,x3,x4,b) 回车输出: a= 1.0e-004* 0 0 0.2083 -0.0000 >>x=0: 50: 12000;y=a (1)+x*a (2)+x.^2*a(3)+x.^3*a(4);plot(x,y) 回车,输出图像: 4.1曾任英特尔公司董事长的摩尔先生早在1965年时,就观察到一件很有趣的现象: 集成电路上可容纳的零件数量,每隔一年半左右就会增长一倍,性能也提升一倍。 因而发表论文,提出了大大有名的摩尔定律(Moore’sLaw),并预测未来这种增长仍会延续下去。 下面数据中,第二行数据为晶片上晶体数目在不同年代与1959年时数目比较的倍数。 这些数据是推出摩尔定律的依据: 年代 1959 1962 1963 1964 1965 增加倍数 1 3 4 5 6 试从表中数据出发,推导线性拟合的函数表达式。 分析: 通过表中数据绘图,可发现趋向于一次线性方程,正好本题要求推到线性拟合函数,直接使用一次的最小二乘法。 程序代码: functionzxec(x,y) y=y'; fori=1: length(x) forj=1: 2 A(i,j)=x(i)^(j-1); end end [L,U]=lu(A'*A); alpha=U\(L\(A'*y)) 在指令窗口输入: x=[1959,1962,1963,1964,1965];y=[1,3,4,5,6];zxec(x,y) 程序输出: alpha= 1.0e+003* -1.6255 0.0008 其中alpha (1)=-1625.5是函数交y轴的焦点值,0.8是一阶函数斜率 5.1用几种不同的方法求积分 的值。 (1)牛顿-莱布尼茨公式; (2)梯形公式;(3)辛卜生公式;(4)复合梯形公式。 程序代码: (1)牛顿-莱布尼茨 由于 =4arctgx|0到1 输入内容: >>4*atan (1) 程序输出: ans= 3.1416 (2)(3)(4) 定义 functions=f(x) s=4/(1+x^2); 输入: a=0;b=1* (2)梯形公式 接着*式输入: (1/2)*(b-a)*(f(b)+f(a)) 输出结果: ans= 3 (3)辛卜生公式 接着*式输入: (1/6)*(b-a)*(f(a)+4*f((a+b)/2)+f(b)) 输出结果: ans= 3.1333 (4)复合梯形公式 接着*式输入: T=0; forn=0: 0.1: 0.9 m=0.1/2*(f(n)+f(n+0.1));T=T+m; end T 输出结果: T= 3.1399 6.1用欧拉公式和四阶龙格-库塔法分别求解下列初值问题; 分析: 此题只需建立欧拉和龙格库塔函数,带入求解即可。 由于题目没有给出步长,我们取步长h=0.1 程序代码: (1)建立名为f.m的文件: functionr=f(t,y) r=(0.9*y)/(1+2*t); Euler方法 functioneuler(a,b,h,alpha) t=a: h: b; n=(b-a)/h; y (1)=alpha; fori=2: n+1 y(i)=y(i-1)+h*f(t(i-1),y(i-1)); sprintf('t=%3.1f,y=%9.7f',t(i),y(i)) end 在命令窗口输入: euler(0,1,0.1,1) 回车得到: ans= t=0.1,y=1.0900000 t=0.2,y=1.1717500 t=0.3,y=1.2470768 t=0.4,y=1.3172249 t=0.5,y=1.3830861 t=0.6,y=1.4453250 t=0.7,y=1.5044519 t=0.8,y=1.5608688 t=0.9,y=1.6148989 t=1.0,y=1.6668064 四阶Runge-Kutta方法 functionrk4(a,b,h,alpha) t=a: h: b; n=(b-a)/h; y (1)=alpha; fori=2: n+1 k1=f(t(i-1),y(i-1)); k2=f(t(i-1)+h/2,y(i-1)+h/2*k1); k3=f(t(i-1)+h/2,y(i-1)+h/2*k2); k4=f(t(i-1)+h,y(i-1)+h*k3); y(i)=y(i-1)+h/6*(k1+2*k2+2*k3+k4); sprintf('t=%3.1f,y=%9.7f',t(i),y(i)) end 在命令窗口输入: rk4(0,1,0.1,1) 回车得到: ans= t=0.1,y=1.0855051 t=0.2,y=1.1634777 t=0.3,y=1.2355334 t=0.4,y=1.3027862 t=0.5,y=1.3660420 t=0.6,y=1.4259056 t=0.7,y=1.4828446 t=0.8,y=1.5372290 t=0.9,y=1.5893579 t=1.0,y=1.6394763 同理: (2)建立名为f.m的文件: functionr=f(t,y) r=-(t*y)/(1+t^2); 在命令窗口分别输入: euler(0,1,0.1,2) 回车 rk4(0,1,0.1,2) 回车 结果如下表 Euler方法 四阶Runge-Kutta方法 t=0.1,y=2.0000000 t=0.2,y=1.9801980 t=0.3,y=1.9421173 t=0.4,y=1.8886645 t=0.5,y=1.8235382 t=0.6,y=1.7505966 t=0.7,y=1.6733644 t=0.8,y=1.5947500 t=0.9,y=1.5169573 t=1.0,y=1.4415285 t=0.1,y=1.9900743 t=0.2,y=1.9611612 t=0.3,y=1.9156523 t=0.4,y=1.8569530 t=0.5,y=1.7888539 t=0.6,y=1.7149854 t=0.7,y=1.6384634 t=0.8,y=1.5617372 t=0.9,y=1.4865879 t=1.0,y=1.4142132 7.1用二分法求下列方程在指定区间[a,b]上的实根近似值: (要求误差不超过0.01) (1)x-sinx-1=0,[a,b]=[1,3]; (2)xsinx=1,[a,b]=[0,1.5]. 分析: 直接将方程带入二分法函数中求解即可。 程序代码: functionx=eff(fun,a,b,eps) ifnargin<4 eps=1e-5; end fa=feval(fun,a);fb=feval(fun,b); iffa*fb>0 disp('[a,b]不是有根区间,请重新调整'); return; end whileabs(b-a)/2>eps x=(a+b)/2;fx=feval(fun,x); iffx*fa<0 b=x;fb=fx; else a=x;fa=fx; end end x=(a+b)/2; 在命令窗口输入: (1)fun=inline('x-sin(x)-1');x=eff(fun,1,3,0.01) 输出结果 x= 1.9297 (2)fun=inline('x*sin(x)-1');x=eff(fun,0,1.5,0.01) 输出结果: x= 1.1191 8.1分别用直接法、雅可比迭代法、赛德尔迭代法求解线性方程组AX=b。 (1) , (2) , 分析: 首先通过diag函数给矩阵A附上初值,然后分别调用雅可比和赛德尔函数即可 程序: (1) M=ones(9,9); N=4*ones(10,10); A=diag(diag(N))+diag(diag(M),1)+diag(diag(M),-1); b=[2;1;1;1;1;1;1;1;1;2]; 1.用直接法输入: formatlong A\b 程序输出: ans= 0.479271991911021 0.082912032355915 0.189********5319 0.160768452982811 0.167846309403438 0.167846309403438 0.160768452982811 0.189********5319 0.082912032355915 0.479271991911021 2.用Jacobi法 程序代码: functionJacobi(A,b,X0,P,error,maxl) [nn]=size(A); X=zeros(n,1); fork=1: maxl forj=1: n X(j)=(b(j)-A(j,[1: j-1,j+1: n])*X0([1: j-1,j+1: n]))/A(j,j); end errX=norm(X-X0,P); X0=X; if(errX disp('迭代次数k,近似解X分别是: ') k X return end end if(errX>=error) disp('请注意: Jacobi迭代次数已经超过最大迭代次数maxl.') end 在命令窗口输入: X0=[0;0;0;0;0;0;0;0;0;0]; Jacobi(A,b,X0,inf,0.001,100); 程序输出: 迭代次数k,精确解X1和近似解X分别是: k= 9 X= 0.47933959960938 .0831******** 0.189********852 0.16108322143555 0.16813659667969 0.16813659667969 0.16108322143555 0.189********852 .0831******** 0.47933959960938 用赛德尔迭代法: 程序代码: functionGS(A,b,X0,P,error,maxl) [nn]=size(A); X=zeros(n,1); fork=1: maxl forj=1: n XX=0; fori=1: n ifi XX=XX+A(j,i)*X(i); end ifi>j XX=XX+A(j,i)*X0(i); end end X(j)=(b(j)-XX)/A(j,j); end errX=norm(X-X0,P); X0=X; if(errX disp('迭代次数k,近似解X分别是: ') k X return end end if(errX>=error) disp('请注意: Gauss-Seidel迭代次数已经超过最大迭代次数maxl.') end 在命令窗口输入: X0=[0;0;0;0;0;0;0;0;0;0]; GS(A,b,X0,inf,0.001,100); 输出结果: k= 6 X= 0.47925853729248 0.08289575576782 0.189********913 0.16064277291298 0.16802297346294 0.16770140314475 0.16085489129182 0.189********114 0.08292683955369 0.47926829011158 (2) U=ones(18,18); V=-2*ones(19,19); W=5*ones(20,20); A=diag(diag(W))+diag(diag(V),1)+diag(diag(U),2)+diag(diag(V),-1)+diag(diag(U),-2); b=[1;zeros(19,1)]; X0=zeros(20,1); 同上获得如下解: 直接法 雅可比 赛德尔 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 0.242121373548085 .0943******** .021********* .0313******** -0.006942557862112 0.004886927715767 0.003586117214438 0.000214823135181 -0.000784526508185 -0.000357861979163 0.000050437638520 0.000106309021306 0.000029232399870 -0.000014343050585 -0.000012667696430 -0.000001464361499 0.000002491258113 0.000001311981447 -0.000000093354238 -0.000000299737984 1.0e+004*-0.28210827857172 1.0e+004*0.48204640630127 1.0e+004*-0.69195375595679 1.0e+004*0.88126694207517 1.0e+004*-1.05291766695873 1.0e+004*1.20130020954019 1.0e+004*-1.32375170735604 1.0e+004*1.41753761486760 1.0e+004*-1.48072923013544 1.0e+004*1.51202411408633 1.0e+004*-1.51082661747071 1.0e+004*1.47723879582980 1.0e+004*-1.41205236153589 1.0e+004*1.31674011324409 1.0e+004*-1.19335996640489 1.0e+004*1.04472588730159 1.0e+004*-0.87353157471502 1.0e+004*0.68530915502268 1.0e+004*-0.47711982856112 1.0e+004*0.27913406015473 0.24229171609600 0.09477250842624 .021********* .0314******** -0.0071474575114
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 课程设计