数值实验报告

1、2008-2009数值实验数值实验报告一、为了逼近飞行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。见下表。1 对这些数据构造三次自然样条插值函数，并画出得到的三次自然样条插值曲线；2 对这些数据构造Lagrange插值多项式，并画出得到的Lagrange插值多项式曲线。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解：程序:【La

2、gr.m】function y=lagr(x0,y0,x)n=length(x0); m=length(x);for i=1:m z=x(i); s=0; for k=1:n p=1; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end【exam1.m】x0=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;y0=1.3 1.5 1.85 2.1 2.6 2

3、.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;x=0.9:0.02:13.3;y1=spline(x0,y0,x);y2=lagr(x0,y0,x);plot(x0,y0,'*',x,y1,'-',x,y2,'r')计算结果：*为所给的数据点，绿色的虚线为三次样条曲线，红色实线为Lagrange插值多项式曲线。二、对于问题 将h=0.025的Euler法，h=0.05的改进的Euler法和h=0.1的4阶经典的Runge-Kutta法在这些方法的公共节点0

4、.1，0.2，0.3，0.4和0.5处进行比较。精确解为：。解：程序：【euler.m】function u=euler(t,h)u(1)=0.5;n=length(t)-1;for i=1:n u(i+1)=u(i)+h*(u(i)-(t(i)2+1);end【euler_gaijin.m】function u=euler_gaijin(t,h)u(1)=0.5;n=length(t)-1;for i=1:n uu(i+1)=u(i)+h*(u(i)-(t(i)2+1); u(i+1)=u(i)+0.5*h*(u(i)-(t(i)2+1)+(uu(i+1)-(t(i+1)2+1);end【r

5、ungekutta.m】function u=rungekutta(t,h)u(1)=0.5;n=length(t)-1;for i=1:20 k1=u(i)-(t(i)2+1; k2=(u(i)+0.5*h*k1)-(t(i)+0.5*h)2+1; k3=(u(i)+0.5*h*k2)-(t(i)+0.5*h)2+1; k4=(u(i)+h*k3)-(t(i)+h)2+1; u(i+1)=u(i)+h/6*(k1+2*k2+2*k3+k4);end【exam2.m】clc;t=0:0.01:2;u=(1+t).2-0.5*exp(t);h1=0.025;t1=0:h1:2;u1=euler(

6、t1,h1);h2=0.05;t2=0:h2:2;u2=euler_gaijin(t2,h2);h3=0.1;t3=0:h3:2;u3=rungekutta(t3,h3);plot(t,u,'g',t1,u1,'-'),pauseplot(t,u,'g',t2,u2,'*-'),pauseplot(t,u,'g',t3,u3,'*-')x1=t(11),u(11); t1(5),u1(5); t2(3),u2(3); t3(2),u3(2);x2=t(21),u(21); t1(9),u1(9);

7、t2(5),u2(5); t3(3),u3(3);x3=t(31),u(31); t1(13),u1(13); t2(7),u2(7); t3(4),u3(4);x4=t(41),u(41); t1(17),u1(17); t2(9),u2(9); t3(5),u3(5);x5=t(51),u(51); t1(21),u1(21); t2(11),u2(11); t3(6),u3(6);计算结果：Euler法：其中的绿色实线为精确解的曲线，蓝色虚线为Euler法计算结果曲线改进Euler法：其中的绿色实线为精确解的曲线，蓝色*线为改进Euler法计算结果曲线Runge-Kutta法：其中的绿色

8、实线为精确解的曲线，蓝色*线为Runge-Kutta法计算结果曲线在公共节点的比较：x1 = 0.10000000000000 0.65741454096218 精确解 0.10000000000000 0.65549823242187 Euler法计算结果 0.10000000000000 0.65730851562500 改进Euler法计算结果 0.10000000000000 0.65741437500000 Runge-Kutta法计算结果x2 = 0.20000000000000 0.82929862091992 精确解 0.20000000000000 0.82533847880

9、729 Euler法计算结果 0.20000000000000 0.82907775662463 改进Euler法计算结果 0.20000000000000 0.82929827599740 Runge-Kutta法计算结果x3 = 0.30000000000000 1.01507059621200 精确解 0.30000000000000 1.00893336727069 Euler法计算结果 0.30000000000000 1.01472539809879 改进Euler法计算结果 0.30000000000000 1.01507005843261 Runge-Kutta法计算结果x4

10、= 0.40000000000000 1.21408765117936 精确解 0.40000000000000 1.20563454915320 Euler法计算结果 0.40000000000000 1.21360789733236 改进Euler法计算结果 0.40000000000000 1.21408690570301 Runge-Kutta法计算结果x5 = 0.50000000000000 1.42563936464994 精确解 0.50000000000000 1.41472636884754 Euler法计算结果 0.50000000000000 1.42501405817

11、677 改进Euler法计算结果 0.50000000000000 1.42563839564822 Runge-Kutta法计算结果三、用Newton迭代法求方程的根时，分别取初始值，； 用Newton迭代法求方程时，分别取初始值，；解： 程序：【newton1.m】clc;clear all;x=zeros(1,20);%x(1)=1.45;x(1)=0.5;for i=1:20 x(i+1)=x(i)-atan(x(i)*(1+x(i)2); if abs(x(i+1)-x(i)<1e-8 end x1=x(i+1);endformat longx1,x'【newton2.

12、m】clc;clear all;x=zeros(1,20);%x(1)=2.0;for i=1:20 x(i+1)=x(i)-(x(i)3-x(i)-3)/(3*x(i)2-1); if abs(x(i+1)-x(i)<1e-8 end x2=x(i+1);endformat longx2,x'【newton1.m】计算结果：（1）时，x1 = NaNans = 1.0e+281 * 0.00000000000000 -0.00000000000000 0.00000000000000 -0.00000000000000 0.00000000000000 -0.000000000

13、00000 0.00000000000000 -0.00000000000000 0.00000000000000 -0.00000000000000 0.00000000000000 -0.00000000000000 7.99254161852540 -Inf NaN NaN NaN NaN NaN NaN NaN迭代发散。（2）时，x1 = 0ans = 0.50000000000000 -0.07955951125101 0.00033530220401 -0.00000000002513 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0迭代到第4步就收敛，解为。【n

14、ewton2.m】计算结果：（1）时，ans = 0 -3.00000000000000 -1.96153846153846 -1.14717596140355 -0.00657937148071 -3.00038907407123 -1.96181817566632 -1.14743022848160 -0.00725624755242 -3.00047318877322 -1.96187864636024 -1.14748519321677 -0.00740250133287 -3.00049244291695 -1.96189248824631 -1.14749777454454 -0.

15、00743597525670 -3.00049690365354 -1.96189569508551 -1.14750068932977 -0.00744373016891迭代发散。（1）时，x2 = 1.67169988165716ans = 2.00000000000000 1.72727272727273 1.67369117369117 1.67170256974750 1.67169988166207 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169

16、988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716 1.67169988165716迭代到第5步就收敛，解为。四、用Jacobi迭代和SOR法分别求解线性方程组（教科书第86页算例2），并验证、输出SOR法的松弛因子w和对应的迭代次数。解： 程序：【data.m】4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0-

17、1 4 -10 -1 4 -10 0 -1 4 -10 0 0 -1 4 -10 0 0 0 -1 4 -10 0 0 0 0 -1 4 -10 0 0 0 0 0 -1 4 -10 0 0 0 0 0 0 -1 4 -10 0 0 0 0 0 0 0 -1 4 -10 0 0 0 0 0 0 0 0 -1 4 -10 0 0 0 0 0 0 0 0 0 -1 4 -10 0 0 0 0 0 0 0 0 0 0 -1 4 -10 0 0 0 0 0 0 0 0 0 0 0 -1 4 -10 0 0 0 0 0 0 0 0 0 0 0 0 -1 4【Jacobi.m】clc;n=15;I=eye

18、(n,n);A=dlmread('data.m');b=zeros(n,1);b(1)=3;b(n)=3;b(2:n-1)=2;D=diag(diag(A);D=inv(D);Bj=I-D*A;f=D*b;err=1;x=zeros(n,1);cnt=0;while abs(err)>1e-6 xx=Bj*x+f; err=norm(xx-x); x=xx; cnt=cnt+1;endformat shortx1=x,cnt【SOR.m】clc;n=15;I=eye(n,n);A=dlmread('bankdata.m');b=zeros(n,1);b(1)=3;b(n)=3