矩阵与数值分析实验报告.docx

《矩阵与数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《矩阵与数值分析实验报告.docx（22页珍藏版）》请在第一文库网上搜索。

1、矩阵与数值分析课程数值实验报告数值实验报告内容：要包含题目、算法公式、完整的程序、正确的数值结果和图形以及相应的误差分析。一、为了逼近飞行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。见下表。对这些数据构造三次自然样条插值函数，并画出得到的三次自然样条插值曲线；2.对这些数据构造Lagrange插值多项式，并画出得到的Lagrange插值多项式曲线。X0.91.31.92.12.63.03.94.44.75.06.0於)1.31.51.852.12.62.72.42.152.052.12.25X7.08.09.210.511.311.612.012.613.013.3於)2.32.25

2、1.951.40.90.70.60.50.40.25解答：1.构造的三次自然样条插值函数;a.算法公式：设给定节点a = x0xl-xn=h及节点上的函数值() = Xz=0,1,构造三次样条插值问题就是构造s(x)S3(2,X,”)，使s()=y=o, 当X*,Z+J时，s(x)表达式为其中4 =%-4(女=0/,-1)，/=4,叫=力。考虑第二边界条件情况下，即s%)/(%)= /；即，边界条件可以表示为2 =22 J01%2 2mn-l+加=3(-1) hn-%-291fn为求解（左=0,1,，使用方程组其中。_3(必一打)%)f”So = ；-Jo423(-i) -i fo -vz J

3、 nw-b.源程序：format long g;xk=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;%节点处的值h=zeros(l,20);%分别设置初始状态为0向量c=zeros(l,20);d=zeros(l,20);g=zeros(l,21);b=zeros(l

4、,21);m2=zeros(21,1);ml=zeros(21,l);%分别设置初始状态为0向量for i=l:1:20h(i)=xk(i+l)-xk(i);% 算出 hkendfor i=2:l:20c(i)=h(i)(h(i)+h(i-1); %计算g 中的两个系数d(i)=h(i-l)/(h(i-l)+h(i);endfor i=2:l:20g(i)=3*(d(i)*(y(i 1 )-y(i)h(i)+c(i)*(y(i)-y(i-1 )h(i-1);endg( 1 )=3*(y(2)-y( 1 )h( 1);%已知的两个边界条件g(21)=3*(y(21)-y(20)/h(20);v=

5、12 2222222222222222222 2J;% 对角元为 2A=diag(v);for i=2:l:20A(i,i+l)=d(i);A(i,i-l)=c(i);b(i)=g(i);endA(l,2)=l;A(21,20)=1;b二g；b(21)=g(21);L,U=lu(A);ml=Lb;m=Um 1;for i=l:l:20;%转为对角阵% LU分解%m即为解出的m+1个系数分段插值x=xk(i):(xk(i+l)-xk(i)/10:xk(i+l);sx(i,:)=(1 -2.*(x-xk(i)(-h(i).*(x-xk(i+ l),(-h(i)(2).*y(i)+ .(l-2.*(

6、x-xk(i+ l).h(i).*(x-xk(i).h(i).(2).*y(i+1)+ .(x-xk(i).*(x-xk(i+1 ).-h(i).(2).*m(i) .(x-xk(i 1 ).*(x-xk ).h ).人(2).*m(i+1);xx(i,:)=x;endsxl=zeros( 1,201);for j= 1:20;forsxl(l,i+l O*(j-l)=sx(j,i)；xxl(l,i+l 0*(j-l)=xx(j,i)；end%画出图形%得到向量mendplot(xxl,sxl);m=mc.运行结果和图形结果：m= 0.5396238492562310.420752301487

7、5381.086802718677961.294941983029580.593399322097993-0.022191145976441-0.503406025873617-0.477075060628503-0.07131619046462220.2623398224869930.08077550666147850.0145581508670924-0.139008110129848-0.335834096469179-0.531829914635186-0.731178228262683-0.492948654289434-0.141335308965742-0.17890047373

8、7137-0.392774881565715-0.553612559217141图形:结果分析三次样条插值不仅有较好的收敛性和稳定性，而且其光滑性也较高。2.对这些数据构造Lagrange插值多项式a.算法公式：设,X,x在,上的+ 1互异点j(x)(x-x0)(x-xy.1)(x-xy+1)(x-x)(xj -Xo)(- -x7-1)(xy -xy+1)(xy-xfl)ij()=-71P (x) = Zj(x)六()b.源程序：首先定义一个lagrang函数，其有一个输出：function y=lagrange(x,y,x)n=length(xO);p=length(x);for i=l:p

9、s=0;z=x(i)；for m=l:nL=l;for d=l :nif m二dL=L*(z-xO(d)/(xO(m)-xO(d);endends=s+L*yO(m);endy(i)=s；end然后，运行如下代码：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=ll.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.25J;p=100

10、;a=0.9;b=13.3;%x20( 1 :n)=a:(b-a)/(n-1 ):b;%d( 1 :p)=a:(b-a)/(p-1 ):b;d=0.9r(b-a)(p-1):13.3;yO=lagrangc(x,y,d);yl=interpl(x,y,d);plot(d,yObd,ylr,)c.结果和图形File Edit View Insert Tools Desktop Window Help5 隰I园 口结果分析对于拉格朗日插值而言，插值节点多导致多项式的次数增多，而高次多项式的震荡次数增多，有可能在非节点处误差变得很大。二、对于问题/二一/+1,or2,Y,w(0) = 0.5将=0.

11、025的Euler法，=0.05的改进的Euler法和z=0.1的4阶经典的RmgeKtta法在这些方法的公共节点0.1, 0.2, 0.3, 0.4和0.5处进行比较。精确解为：(。=(1 + /)2().5-。解答：a.算法公式：将节点取为b ctk = a + kh, h =, Z = 0,l,2,UM =uk + /(，,0)力,jtk(/0)= 0Euler 法um = w + hf(tk.uk)w =uk+hf (，4)，改进的Euler法h-%=k +j(L% ) + (+1，以+】)，乙w(r0) = w04阶经典的Runge-Kutta法hk+ = % + ( + 242 +

12、 2% + k4),6k= fW,k),1,hh 1、,k? =+展4 +/),jrzhh 1 %3 = /(4 +不，& +5”2)，k2 = f(tk +h.uk + 秘3)，.b.源程序：format long g;a=0; b=2; hl=0.025;h2=0.05;h3=0.1; h=0.1;% 己知的初始值for k=l:l:21t(k)=a+(k-l)*h3;true(k)=-0.5*exp(t(k)+( 1 +t(k)A2;% 6 points 根据解出的真值计算endyl =0.5;tl(l)=0;for k=l :1:2000/25tl(k)=a+(k-l)*hl;yl(k

13、+l)=y l(k)+h 1 *(y l(k)-tl(k)*tl(k)+1);endtl(81)=2;y2 =0.5;t2(l)=0;for k=l:l:40t2(k+1 )=a+(k-0)*h2;y(k+1 )=y2(k)+h2*(y2(k)-t2(k)*t2(k)+1);y2(k+ I)=y2(k)+.5*h2*(y2(k)-t2(k)*t2(k)+ l)(y(k+1 )-t2(k+ l)*t2(k+1)+1);endy3(l)=0.5;% 经典 Runge-Kutta 法%t3(l)=0;for i=l:l:20t3(i)=a+(i-l)*h3;kl=y3(i)-(t3(i)*(t3(i)+l;k2=(y3(i)+0.5*h3*k 1 )-(t3(i)+O.5*h3)*(t3(i)+O.5*h3)+1;k3=(y3 +0.5*h3*k2)(t3 +O.5*h3)*(t3 +0.5*h3)+1;k4=(y3 +h3*k3)(t3 +h3)*(t3 +h3)+1;y3(i+l)=y3 +h3/6*(kl+2*k2+2*k3+k4);endt3(21)=2;y3;subplot(2,2,l)plot(t,true,