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

1、矩阵与数值分析实验报告学院：土木市政工程 姓名：徐博书 学号：教师：张宏伟 班级：2班一、为了逼近飞行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。见下表。.对这些数据构造三次自然样条插值函数，并画出得到的三次自然样条插值曲线；.对这些数据构造 Lagrang插值多项式，并画出得到的 Lagrang插值多项式曲线。xf(x)xf(x)解：.三次样条插值函数程序代码及运行结果clear all ;clc;x=；y=;n=length(y);h=zeros(1,n-1);for i=1:n-1h(i)=x(i+1)-x(i);endr=zeros(1,n-2);r(i)=h(i)/(h(i

2、)+h(i-1);endu=zeros(1,n-2);for i=2:n-1u(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-1g(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-1for j=1:nif i=jA(i,j)=2;A(i,j-1)=r(i

3、);A(i,j+1)=u(i);endendendm=Ag;for i=1:n-1z=x(i):x(i+1);s=(h(i)+2.*(z-x(i)/h(i)A3).*(z-x(i+1)A2).*y(i)+(h(i)-2.*(z-x(i+1)/h(i)A3).*(z-x(i)A2).*y(i+1)+(z-x(i).*(z-x(i+1)A2).*m(i)/(h(i)A2)+(z-x(i+1).*(z-x(i)A2).*m(i+1)/(h(i)A2);plot(x(i),y(i),kp ,z,s);hold onendgrid on;程序运行结果：M=,；. Lagrange插值多项式程序代码及运行

4、结果clear all ;clc;x=;y=;n=length(y);l=ones(1,n);for i=1:nfor j=1:nif i=jl(i)=l(i);else l(i)=l(i)/(x(i)-x(j);endendendl=l.*y;for i=1:n-1p=zeros();z=x(i):x(i+1);for j=1:ns=ones();for q=1:nif j=qs=s.*(z-x(q);endendp=p+l(j).*s;endplot(x(i),y(i),ko ,z,p);hold onendgrid on;程序运行结果：6二、对于问题 TOC o 1-5 h z 22 H

5、YPERLINK l bookmark4 o Current Document u ut1,0t2,、u00.5将卜=的Euler法，卜=的改进的Euler法和=的4阶经典的Runge-Kutta法在这些方法的公共节点，和处进行比较。精确解为： u t 1 t 2 0.5et o1.卜=的Euler法完整MATLA翼序新建M-文件，建立euler1函数，保存至系统默认路径。function x,y=euler1(dyfun,xspan,y0,h)x=xspan:h:xspan(2);y(1)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n

6、),y(n);endx=x;y=y;在MATLA瑜令窗口中输入dyfun=inline(u-tA2+1);t,u=euler1(dyfun,0,2,;t,u回车运行得到ans =的改进的Euler 法完整MATLA叁序新建M-文件，建立euler2函数，保存至系统默认路径。function x,y=euler2(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=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)

7、=y(n)+h*(k1+k2)/2;endx=x;y=y;在MATLA瑜令窗口中输入dyfun=inline(u-tA2+1);t,u=euler2(dyfun,0,2,;t,u回车运行得到ans =的 4 阶经典的 Runge-Kutta 法完整MATLA叁序新建M-文件，建立rungekutta函数，保存至系统默认路径。function x,y=rungekutta(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=feval(dyfun,x(n),y(n);k2=feval(dyfun,x(n)+h/

8、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;在MATLA瑜令窗口中输入dyfun=inline(u-tA2+1);t,u=rungekutta(dyfun,0,2,;t,u回车运行得到ans =4. 精确解在MATLA瑜令窗口中输入t=0:2;u=(1+t).A.*exp(t);t,u回车运行得到ans =5.求常微分方程的方法对比计算值对比h-白Euler 法h=的改进的Eule

9、r 法h= 的 4 阶Runge-Kutta 法精确解本例说明卜二的4阶Runge-Kutta法与真值最为接近。三、用Newton迭代法求方程f x arctan x 0的根时，分别取初始值 x0 1.45 , x0 0.5 ;用Newton迭代法求方程fx x3 x 3 0时，分别取初始值x0 0.0, x0 2.0；xA3-x-33*xA2-1解：.求解 f x arctan x 0取初始值完整MATLA叁序新建M-文件，建立newton函数，保存至系统默认路径。function nx,k=newton(fname,dfname,x0,ep,nmax)if nargin5 nmax=500

10、; endif narginep&k In newton at 9nx = k =500500 次。说明取初值0 500 次。取初始值2在MATLA瑜令窗口中输入fname=inline( xA3-x-3);dfname=inline(3*xA2-1);nx,k=newton(fname,dfname,2)回车运行可得nx =4此题说明了初值的选取对newton 迭代法很重要，选取适当的初值会很快收敛到满足要求的解，选取不当的初值有可能不收敛或收敛很慢。四、用Jacobi迭代和SORt分别求解线性方程组(教科书第86页算例2),并验证、输出SOR去的松弓因子w和对应的迭代次数。解：1、 Jac

11、obi 迭代程序：%第四题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-1for j=1:nif i=jA(i,j)=4;A(i,j+1)=-1;A(i,j-1)=-1;endendend%建立矩阵bb=zeros(1,n);b(1)=3; b(n)=3;for i=2:n-1b(i)=2;end%inv(d)*(L+U)d=(1/4)*eye(n);LU=-1*A+4*eye(n);D=d*LU;%inv(d)*bB=d*(b);%

12、开始迭代q=0;x=zeros(n,1);e=1;while e10A(-6)x=D*x+B;e=norm(x-ones(1,n);q=q+1;if q=40break ;endenddisp( 迭代次数为 );disp( 迭代结果为 );x程序运行结果：当输入A的阶次为4时，迭代次数为17,迭代结果为x=1 .00002、SORt迭代程序：蹴四题SORI迭代clear all ;clc;n=input(系数矩阵 A勺阶数n=);w=input( 松驰因子 w=);D=4*eye(n);L=zeros(n,n);for i=2:nfor j=1:nif i=jL(i,j-1)=1;endendendU=zeros(n,n);for i=1:n-1for j=1:nif i=jU(i,j+1)=1;endendb=zeros(1,n);b(1)=3; b(n)=3;for i=2:n-1b(i)=2;endL1=inv(D-