1、插值法1. 下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间0,64上作图.(1)用这 9 个点作 8 次多项式插值 Ls(x).(2)用三次样条(第一边界条件)程序求 S(x).从得到结果看在0,64上，哪个插值更精确；在区间0,1上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=0 1 4 9 16 25 36 49 64;y1=0 1 2 3 4 5 6 7 8;n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i);for k=1
2、:i-1l=l*(x-x1(k)/(x1(i)-x1(k);endfor k=i+1:nl=l*(x-x1(k)/(x1(i)-x1(k);endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式 Ls(x).输出结果为Ls =-24221063/63504000*x2+95549/72072*x-1/3048192000*x8-2168879/435456000*x4+19/283046400*x7+657859/10886400*x3+33983/152409600*x5-13003/2395008000*x6(2)三次样条插值，程序如下x1=0 1 4 9 16 2
3、5 36 49 64;y1=0 1 2 3 4 5 6 7 8;x2=0:1:64;y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数S=p(1)+p(2)*x+p(3)*x2+p(4)*x3 %得到 S(x)输出结果为：S = 2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x2+999337332656867/1125899906842624*x3(3)在区间0,64
4、上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),b,x2,y2,r,x2,y3,y)蓝色曲线为 y= 函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线0 10 20 30 40 50 60 70-20020406080100可以看到蓝色曲线与黄色曲线几乎重合，因此在区间0,64上三次样条插值更精确。在0,1区间上由上图看不出差别，不妨代入几组数据进行比较 ，取 x4=0:0.2:1x4=0:0.2:1;sqrt(x4) %准确值subs(Ls,x,x4) %拉格朗日插值spline(x1,y1,x4) %三次样条插值运行结果为ans =0 0.4472
5、 0.6325 0.7746 0.8944 1.0000ans =0 0.2504 0.4730 0.6706 0.8455 1.0000ans =0 0.2429 0.4630 0.6617 0.8403 1.0000从这几组数值上可以看出在0,1区间上，拉格朗日插值更精确。数据拟合和最佳平方逼近2. 有实验给出数据表x 0.0 0.1 0.2 0.3 0.5 0.8 1.0y 1.0 0.41 0.50 0.61 0.91 2.02 2.46试求 3 次、4 次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。解：(1)三次拟合，程序如下s
6、ym x;x0=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y0=1.0 0.41 0.50 0.61 0.91 2.02 2.46;cc=polyfit(x0,y0,3);S3=cc(1)+cc(2)*x+cc(3)*x2+cc(4)*x3 %三次拟合多项式xx=x0(1):0.1:x0(length(x0);yy=polyval(cc,xx);plot(xx,yy,-);hold on;plot(x0,y0,x);xlabel(x);ylabel(y);运行结果S3 =-7455778416425083/1125899906842624+1803512222945437/140
7、737488355328*x-655705280524945/140737488355328*x2+4172976892199509/4503599627370496*x3图像如下0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.511.522.5xy(2)4 次多项式拟合sym x;x0=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y0=1.0 0.41 0.50 0.61 0.91 2.02 2.46;cc=polyfit(x0,y0,4);S3=cc(1)+cc(2)*x+cc(3)*x2+cc(4)*x3+cc(5)*x4xx=x0(1):
8、0.1:x0(length(x0);yy=polyval(cc,xx);plot(xx,yy,r);hold on;plot(x0,y0,x);xlabel(x);ylabel(y);运行结果S3 = 3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931990070637/281474976710656*x2-5965836931688425/1125899906842624*x3+8491275464650307/9007199254740992*x4图像如下0 0.1 0.2 0.3 0.4
9、 0.5 0.6 0.7 0.8 0.9 100.511.522.5xy(3)另一个拟合曲线,新建一个 M-file程序如下：function C,L=lagran(x,y)w=length(x);n=w-1;L=zeros(w,w);for k=1:n+1V=1;for j=1:n+1if k=jV=conv(V,poly(x(j)/(x(k)-x(j);endendL(k,:)=V;endC=y*L在命令窗口中输入以下的命令：x=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y=1.0 0.41 0.50 0.61 0.91 2.02 2.46;cc=polyfit(x,y,4)
10、;xx=x(1):0.1:x(length(x);yy=polyval(cc,xx);plot(xx,yy,r);hold on;plot(x,y,x);xlabel(x);ylabel(y);x=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y=0.94 0.58 0.47 0.52 1.00 2.00 2.46; %y 中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据C,L=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);hold on;plot(xx,yy,b,x,y,.);图像如下0 0.1 0.2 0.3 0.4 0.5 0.
11、6 0.7 0.8 0.9 100.511.522.5xy解线性方程组的直接解法3. 线性方程组 Ax=b 的 A 及 b 为A= ，b= ，则解 x= .用 MATLAB 内部函数求 detA10 7 8 77 5 6 58 6 10 97 5 9 10 32233331 1111及 A 的所有特征值和 cond(A)2.若令A+A= ，求解(A+A)(x+ x)=b,输出向量 x 和10 7 8.1 7.27.08 5.04 6 58 5.98 9.89 96.99 5 9 9.98|x|2.从理论结果和实际计算两方面分析线性方程组 Ax=b 解得相对误差|x|2/|x|2 及 A 的相对
12、误差|A|2/|A|2 的关系.解：(1)程序如下clear;A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;det(A)cond(A,2)eig(A)输出结果为ans =1ans =2.9841e+003ans =0.01020.84313.858130.2887(2)程序如下A=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98;b=32 23 33 31;x0=1 1 1 1;x=Ab %扰动后方程组的解x1=x-x0 %x 的值norm(x1,2) %x 的 2-范数运行结果为x =-9.586318.
13、3741-3.22583.5240x1 =-10.586317.3741-4.22582.5240ans =20.9322程序如下A=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98;A0=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;b=32 23 33 31;x0=1 1 1 1;x=Ab;x1=x-x0;norm(x1,2);A1=A-A0 ; %A 的值norm(x1,2)/norm(x0,2) % |x|2/|x|2 的值norm(A1,2)/norm(A0,2) %|A|2/|A|2 的值输出结果为
14、ans =10.4661ans =0.0076可见 A 相对误差只为 0.0076，而得到的结果 x 的相对误差就达到了10.4661，该方程组是病态的，A 的条件数为 2984.1 远远大于 1，当 A 只有很小的误差就会给结果带来很大的影响。非线性方程数值解法4. 求下列方程的实根(1) x2-3x+2-ex=0;(2) x3+2x2+10x-20=0.要求：(1)设计一种不动点迭代法，要使迭代序列收敛，然后再用斯特芬森加速迭代，计算到|x(k)-x(k-1)|1E10break;endif abs(x-x0)Epsbreak;endendi,x运行结果：f =Inline functio
15、n:f(x) = (-2*x2-10*x+20)1/3x=0.166666666666676.09259259259259-38.38843164151806-8.478196837919431e+002-4.763660785374071e+005-1.512815059604763e+011i =6x =-1.512815059604763e+011迭代 6 次后 x 的值大得令人吃惊,表明构造的式子并不收敛.也无法构造出收敛的不动点公式牛顿迭代法，程序如下：format long;x=sym(x);f=sym(x3+2*x2+10*x-20);df=diff(f,x);FX=x-f/df
16、;Fx=inline(FX);disp(x=);x1=0.5;disp(x1);Eps=1E-8;i=0;while 1x0=x1;i=i+1;x1=feval(Fx,x1);disp(x1);if abs(x1-x0)Epsbreak;endendi,x1运行结果：x=1.500000000000001.373626373626371.368814819623961.368808107834411.36880810782137i =4x1 =1.36880810782137比较三种方法，牛顿法的收敛性比较好，相比不动点迭代法要构造出收敛的公式比较难，牛顿法迭代次数也较少，收敛速度快，只是对初
17、值的要求很高，几种方法各有利弊，具体采用哪种也需因题而异。常微分方程初值问题数值解法5. 给定初值问题y=-50y+50x2+2x, ;01y(0)=1/3;用经典的四阶 R-K 方法解该问题，步长分别取 h=0.1,0.025,0.01,计算并打印 x=0.1i(i=0,1,10)各点的值，与准确值 y(x)=1/3e(-50x)+x2 比较。解：取步长 h=0.1，程序如下：%经典的四阶 R-K 方法clear;F=-50*y+50*x2+2*x;a=0;b=1;h=0.1;n=(b-a)/0.1;X=a:0.1:b;Y=zeros(1,n+1);Y(1)=1/3;for i=1:nx=X
18、(i);y=Y(i);K1=h*eval(F);x=x+h/2;y=y+K1/2;K2=h*eval(F);y=Y(i)+K2/2;K3=h*eval(F);x=X(i)+h;y=Y(i)+K3;K4=h*eval(F);Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end%准确值temp=;f=dsolve(Dy=-50*y+50*x2+2*x,y(0)=1/3,x);df=zeros(1,n+1);for i=1:n+1temp=subs(f,x,X(i);df(i)=double(vpa(temp);enddisp( 步长 四阶经典 R-K 法 准确值);disp(X,
19、Y,df);运行结果：步长 四阶经典 R-K 法 准确值1.0e+010 *0 0.00000000003333 0.000000000033330.00000000001000 0.00000000046055 0.000000000001220.00000000002000 0.00000000630625 0.000000000004000.00000000003000 0.00000008640494 0.000000000009000.00000000004000 0.00000118436300 0.000000000016000.00000000005000 0.00001623
20、545110 0.000000000025000.00000000006000 0.00022256067134 0.000000000036000.00000000007000 0.00305093542778 0.000000000049000.00000000008000 0.04182323921740 0.000000000064000.00000000009000 0.57332690347809 0.000000000081000.00000000010000 7.85935630083771 0.00000000010000%画图观察结果figure;plot(X,df,k*,
21、X,Y,-r);grid on;title(四阶经典 R-K 法解常微分方程);legend(准确值,四阶经典 R-K 法);0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1012345678x 1010 值值值值R-K值值值值值值值值值值值值值值R-K值当 x 值接近 1 的时候，偏离准确值太大。当步长 h=0.025 时，将上面程序中的 h 改为 0.025 即可，运行结果：步长 四阶经典 R-K 法 准确值0 0.33333333333333 0.333333333333330.10000000000000 0.10313524034288 0.012245
22、982333030.20000000000000 0.04428527327599 0.040015133309920.30000000000000 0.05196795755507 0.090000101967440.40000000000000 0.09395731149439 0.160000000687050.50000000000000 0.16034531435757 0.250000000004630.60000000000000 0.24808570130557 0.360000000000030.70000000000000 0.35624188472758 0.490000
23、000000000.80000000000000 0.48452590661627 0.640000000000000.90000000000000 0.63284923300751 0.810000000000001.00000000000000 0.80118464374206 1.00000000000000图像如下：0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.50.60.70.80.91 值 值 值 值 R-K值 值 值 值 值 值 值值 值 值值 值 值 值 R-K值当步长 h=0.01 时，将上面程序中的 h 改为 0
24、.01 即可，运行结果：步长 四阶经典 R-K 法 准确值0 0.33333333333333 0.333333333333330.10000000000000 0.20235720486111 0.012245982333030.20000000000000 0.12881700190791 0.040015133309920.30000000000000 0.09799182667850 0.090000101967440.40000000000000 0.10094946775024 0.160000000687050.50000000000000 0.13227011975470 0.
25、250000000004630.60000000000000 0.18866520287199 0.360000000000030.70000000000000 0.26813930278431 0.490000000000000.80000000000000 0.36948166028319 0.640000000000000.90000000000000 0.49195762199475 0.810000000000001.00000000000000 0.63512142167910 1.000000000000000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.50.60.70.80.91 值 值 值 值 R-K值 值 值 值 值 值 值值 值 值值 值 值 值 R-K值
