《矩阵与数值分析》实验报告.doc

矩阵与数值分析实验报告学院：土木市政工程 姓名：徐博书 学号：20806075 教师：张宏伟 班级：2班一、为了逼近飞行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。见下表。1 对这些数据构造三次自然样条插值函数，并画出得到的三次自然样条插值曲线；2 对这些数据构造Lagrang插值多项式，并画出得到的Lagrang插值多项式曲线。x0.91.31.92.12.63.03.94.44.75.06.0f(x)1.31.51.852.12.62.72.42.152.052.12.25 x7.08.09.210.511.311.612.012.613.013.3 f(x)2.32.251.951.40.90.70.60.50.40.25解:1.三次样条插值函数程序代码及运行结果clear all;clc;x=0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3;y=1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25;n=length(y);h=zeros(1,n-1);for i=1:n-1 h(i)=x(i+1)-x(i);endr=zeros(1,n-2);for i=2:n-1 r(i)=h(i)/(h(i)+h(i-1);endu=zeros(1,n-2);for i=2:n-1 u(i)=h(i-1)/(h(i)+h(i-1);endg=zeros(1,n);g(1)=3*(y(2)-y(1)/h(1);g(n)=3*(y(n)-y(n-1)/h(n-1);for i=2:n-1 g(i)=3*(u(i)*(y(i-1)-y(i)/h(i)+r(i)*(y(i)-y(i-1)/h(i-1);endA=zeros(n,n);A(1,1)=2; A(1,2)=1; A(n,n)=2; A(n,n-1)=1;for i=2:n-1 for j=1:n if i=j A(i,j)=2; A(i,j-1)=r(i); A(i,j+1)=u(i); end endendm=Ag;for i=1:n-1 z=x(i):0.01:x(i+1);s=(h(i)+2.*(z-x(i)/h(i).3).*(z-x(i+1).2).*y(i)+(h(i)-2.*(z-x(i+1)/h(i).3).*(z-x(i).2).*y(i+1)+(z-x(i).*(z-x(i+1).2).*m(i)/(h(i).2)+(z-x(i+1).*(z-x(i).2).*m(i+1)/(h(i).2); plot(x(i),y(i),kp,z,s); hold on endgrid on;程序运行结果：M= 0.4217,0.6566,-2.6655,2.222,-1.0192,0.5411,0.1305,0.5379,-0.1994,0.2599,-0.0709,0.0238,-0.0243,0.0234,-0.1913,1.8850,-0.8150,0.1235, -0.0199,0.4008,-0.9504； 2. Lagrange插值多项式程序代码及运行结果 clear all;clc;x=0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3;y=1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25;n=length(y);l=ones(1,n);for i=1:n for j=1:n if i=j l(i)=l(i); else l(i)=l(i)/(x(i)-x(j); end endendl=l.*y;for i=1:n-1 p=zeros(); z=x(i):0.01:x(i+1); for j=1:n s=ones(); for q=1:n if j=q s=s.*(z-x(q); end end p=p+l(j).*s; end plot(x(i),y(i),ko,z,p); hold onendgrid on;程序运行结果：二、对于问题 将h=0.025的Euler法，h=0.05的改进的Euler法和h=0.1的4阶经典的Runge-Kutta法在这些方法的公共节点0.1，0.2，0.3，0.4和0.5处进行比较。精确解为：。1. h=0.025的Euler法 完整MATLAB程序新建M-文件，建立euler1函数，保存至系统默认路径。function x,y=euler1(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n);endx=x;y=y;在MATLAB命令窗口中输入dyfun=inline(u-t2+1);t,u=euler1(dyfun,0,2,0.5,0.025);t,u回车运行得到ans = 0 0.5000 0.0250 0.5375 0.0500 0.5759 0.0750 0.6153 0.1000 0.6555 0.1250 0.6966 0.1500 0.7387 0.1750 0.7816 0.2000 0.8253 0.2250 0.8700 0.2500 0.9155 0.2750 0.9618 0.3000 1.0089 0.3250 1.0569 0.3500 1.1057 0.3750 1.1553 0.4000 1.2056 0.4250 1.2568 0.4500 1.3087 0.4750 1.3613 0.5000 1.4147 0.5250 1.4688 0.5500 1.5237 0.5750 1.5792 0.6000 1.6354 0.6250 1.6923 0.6500 1.7498 0.6750 1.8080 0.7000 1.8668 0.7250 1.9263 0.7500 1.9863 0.7750 2.0469 0.8000 2.1080 0.8250 2.1697 0.8500 2.2320 0.8750 2.2947 0.9000 2.3579 0.9250 2.4216 0.9500 2.4858 0.9750 2.5503 1.0000 2.6153 1.0250 2.6807 1.0500 2.7465 1.0750 2.8126 1.1000 2.8790 1.1250 2.9457 1.1500 3.0127 1.1750 3.0800 1.2000 3.1475 1.2250 3.2152 1.2500 3.2830 1.2750 3.3510 1.3000 3.4192 1.3250 3.4874 1.3500 3.5557 1.3750 3.6240 1.4000 3.6924 1.4250 3.7607 1.4500 3.8289 1.4750 3.8971 1.5000 3.9651 1.5250 4.0330 1.5500 4.1007 1.5750 4.1681 1.6000 4.2353 1.6250 4.3022 1.6500 4.3687 1.6750 4.4349 1.7000 4.5006 1.7250 4.5659 1.7500 4.6306 1.7750 4.6948 1.8000 4.7585 1.8250 4.8214 1.8500 4.8837 1.8750 4.9452 1.9000 5.0060 1.9250 5.0659 1.9500 5.1249 1.9750 5.1829 2.0000 5.24002.h=0.05的改进的Euler法完整MATLAB程序新建M-文件，建立euler2函数，保存至系统默认路径。function x,y=euler2(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1); y(n+1)=y(n)+h*(k1+k2)/2;endx=x;y=y;在MATLAB命令窗口中输入dyfun=inline(u-t2+1);t,u=euler2(dyfun,0,2,0.5,0.025);t,u回车运行得到ans = 0 0.5000 0.0500 0.5768 0.1000 0.6573 0.1500 0.7414 0.2000 0.8291 0.2500 0.9202 0.3000 1.0147 0.3500 1.1126 0.4000 1.2136 0.4500 1.3178 0.5000 1.4250 0.5500 1.5352 0.6000 1.6482 0.6500 1.7639 0.7000 1.8822 0.7500 2.0030 0.8000 2.1261 0.8500 2.2514 0.9000 2.3789 0.9500 2.5082 1.0000 2.6393 1.0500 2.7720 1.1000 2.9061 1.1500 3.0415 1.2000 3.1779 1.2500 3.3152 1.3000 3.4531 1.3500 3.5913 1.4000 3.7298 1.4500 3.8682 1.5000 4.0063 1.5500 4.1437 1.6000 4.2803 1.6500 4.4156 1.7000 4.5494 1.7500 4.6814 1.8000 4.8112 1.8500 4.9384 1.9000 5.0627 1.9500 5.1836 2.0000 5.30073.h=0.1的4阶经典的Runge-Kutta法完整MATLAB程序新建M-文件，建立rungekutta函数，保存至系统默认路径。function x,y=rungekutta(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;Endx=x;y=y;在MATLAB命令窗口中输入dyfun=inline(u-t2+1);t,u=rungekutta(dyfun,0,2,0.5,0.025);t,u回车运行得到ans = 0 0.5000 0.1000 0.6574 0.2000 0.8293 0.3000 1.0151 0.4000 1.2141 0.5000 1.4256 0.6000 1.6489 0.7000 1.8831 0.8000 2.1272 0.9000 2.3802 1.0000 2.6409 1.1000 2.9079 1.2000 3.1799 1.3000 3.4553 1.4000 3.7324 1.5000 4.0092 1.6000 4.2835 1.7000 4.5530 1.8000 4.8152 1.9000 5.0670 2.0000 5.30554. 精确解在MATLAB命令窗口中输入t=0:0.1:2;u=(1+t).2-0.5.*exp(t);t,u回车运行得到ans = 0 0.5000 0.1000 0.6574 0.2000 0.8293 0.3000 1.0151 0.4000 1.2141 0.5000 1.4256 0.6000 1.6489 0.7000 1.8831 0.8000 2.1272 0.9000 2.3802 1.0000 2.6409 1.1000 2.9079 1.2000 3.1799 1.3000 3.4554 1.4000 3.7324 1.5000 4.0092 1.6000 4.2835 1.7000 4.5530 1.8000 4.8152 1.9000 5.0671 2.0000 5.30555. 求常微分方程的方法对比计算值对比h=0.025的Euler法h=0.05的改进的Euler法h=0.1的4阶Runge-Kutta法精确解0.10.65550.65730.65740.65740.20.82530.82910.82930.82930.31.00891.01471.01511.01510.41.20561.21361.21411.21410.51.41471.4250 1.42561.4256本例说明h=0.1的4阶Runge-Kutta法与真值最为接近。三、用Newton迭代法求方程的根时，分别取初始值，； 用Newton迭代法求方程时，分别取初始值，； x3-x-3 3*x2-1解：1.求解取初始值1.45完整MATLAB程序新建M-文件，建立newton函数，保存至系统默认路径。function nx,k=newton(fname,dfname,x0,ep,nmax)if nargin5 nmax=500;endif narginep&k In newton at 9nx = -0.0074k = 500说明取初值0时已达到迭代的上限次数500次。取初始值2在MATLAB命令窗口中输入 fname=inline( x3-x-3);dfname=inline(3*x2-1);nx,k=newton(fname,dfname,2)回车运行可得 nx = 1.6717 k = 4此题说明了初值的选取对newton迭代法很重要，选取适当的初值会很快收敛到满足要求的解，选取不当的初值有可能不收敛或收敛很慢。4、 用Jacobi迭代和SOR法分别求解线性方程组（教科书第86页算例2），并验证、输出SOR法的松弛因子w和对应的迭代次数。解： 1、Jacobi迭代程序：%第四题Jacobi迭代n=input(系数矩阵A的阶数n=);%建立系数矩阵AA=zeros(n,n);A(1,1)=4; A(1,2)=-1; A(n,n)=4; A(n,n-1)=-1;for i=2:n-1 for j=1:n if i=j A(i,j)=4; A(i,j+1)=-1; A(i,j-1)=-1; end endend%建立矩阵bb=zeros(1,n);b(1)=3; b(n)=3;for i=2:n-1 b(i)=2;end%inv(d)*(L+U)d=(1/4)*eye(n);LU=-1*