插值法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: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*x^2+95549/72072*x-1/3048192000*x^8-2168879/435456000 *x^4+19/283046400*x^7+657859/10886400*x^3+33983/152409600*x^5-13003/2395008 000*x^6(2)三次样条插值，程序如下x1=[0 1 4 9 16 25 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)*x^2+p(4)*x^3 %得到S(x) 输出结果为： S =23491/304472833/8*x+76713/*x^2+6867/42624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线010203040506070-2020406080100可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。

在[0,1]区间上由上图看不出差别，不妨代入几组数据进行比较 ，取x4=[0:0.2:1]x4=[0:0.2:1];sqrt(x4) %准确值subs(Ls,'x',x4) %拉格朗日插值spline(x1,y1,x4) %三次样条插值运行结果为ans =0 0.4472 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)三次拟合，程序如下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,3);S3=cc(1)+cc(2)*x+cc(3)*x^2+cc(4)*x^3 %三次拟合多项式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 =-25083/42624+45437/5328*x-4945/5328*x^2+99509/70496*x^3图像如下y00.10.20.30.40.50.60.70.80.91x(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)*x^2+cc(4)*x^3+cc(5)*x^4xx=x0(1):0.1:x0(length(x0));yy=polyval(cc,xx);plot(xx,yy,'r');hold on;plot(x0,y0,'x');xlabel('x');ylabel('y');运行结果S3 =96215//0656*x+70637/0656*x^2-88425/42624*x^3+50307/40992*x^4图像如下y00.10.20.30.40.50.60.70.80.91x(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);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,'.');图像如下y00.10.20.30.40.50.60.70.80.91x解线性方程组的直接解法3.线性方程组Ax=b的A及b为A=，b=，则解x=.用MATLAB部函数求det A及A 的所有特征值和cond(A)2.若令A+δA=，求解(A+δA)(x+δx)=b,输出向量x和||δx||2.从理论结果和实际计算两方面分析线性方程组Ax=b解得相对误差||δx||2/||x||2及A的相对误差||δ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=A\b %扰动后方程组的解x1=x-x0 %δx的值norm(x1,2) %δx的2-数运行结果为x =-9.586318.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=A\b;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的值输出结果为ans =10.4661ans =0.0076可见A相对误差只为0.0076，而得到的结果x的相对误差就达到了10.4661，该方程组是病态的，A的条件数为2984.1远远大于1，当A只有很小的误差就会给结果带来很大的影响。

非线性方程数值解法4.求下列方程的实根(1)x^2-3x+2-e^x=0;(2)x^3+2x^2+10x-20=0.要求：(1)设计一种不动点迭代法，要使迭代序列收敛，然后再用斯特芬森加速迭代，计算到|x(k)-x(k-1)|<10^(-8)为止。

(2)用牛顿迭代，同样计算到|x(k)-x(k-1)|<10^(-8)。

输出迭代初值及各次迭代值和迭代次数k ，比较方法的优劣。

解:(1)先用画图的方法估计根的围ezplot('x^2-3*x+2-exp(x)'); grid on;-6-4-20246-200-150-100-5050xx 2-3 x+2-exp(x)可以估计到方程的根在区间（0,1）;选取迭代初值为x0=0.5； 构造不动点迭代公式x(k+1)=( x(k)^2+2-e^x(k))/3; ψ(x)= ( x^2+2-e^x)/3;当0<x<1时，0<ψ(x)<1; ψ’(x)<`1;因此迭代公式收敛。

程序如下：format long;f=inline('(x^2+2-exp(x))/3')disp('x=');x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while 1x0=x;i=i+1;x=feval(f,x);disp(x);if abs(x-x0)<Epsbreak;endendi,x运行结果为f =Inline function:f(x) = (x^2+2-exp(x))/3x=0.9960.8370.4130.4940.5090.2190.4830.5250.7960.8330.0460.5640.7670.501i =14x =0.501斯特芬森加速法，程序如下:format long;f=inline('x-((x^2+2-exp(x))/3-x)^2/((((x^2+2-exp(x))/3)^2+2-exp((x^2+2-exp(x))/3))/3-2*(x^2+2-exp(x))/3+x)');disp('x=');x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while 1x0=x;i=i+1;x=feval(f,x);disp(x);if abs(x-x0)<Epsbreak;endendi,x运行结果为x=0.5790.9810.9860.986i =4x =0.986牛顿迭代法，程序如下：format long;x=sym('x');f=sym('x^2-3*x+2-exp(x)');df=diff(f,x);FX=x-f/df;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=0.0000.8290.4710.9680.986i =4x1 =0.986(2) 先用画图的方法估计根的围ezplot('x^3+2*x^2+10*x-20'); grid on;-6-4-20246-200-100100200300400xx 3+2 x 2+10 x-20根大约在区间（1,2）；选取初值x0=1.5;构造不动点迭代公式x(k+1)=(-2x(k)^2-10x(k)+20)^1/3; ψ(x)=(-2x^2-10x+20)^1/3; 程序如下：format long;f=inline('(-2*x^2-10*x+20)^1/3') disp('x='); x=feval(f,1.5); disp(x) Eps=1E-8; i=1; while 1 x0=x; i=i+1;x=feval(f,x); disp(x);if abs(x-x0)>1E10 break;if abs(x-x0)<Epsbreak;endendi,x运行结果：f =Inline function:f(x) = (-2*x^2-10*x+20)^1/3x=0.6676.259-38.806-8.9431e+002-4.4071e+005-1.4763e+011i =6x =-1.4763e+011迭代6次后x的值大得令人吃惊,表明构造的式子并不收敛. 也无法构造出收敛的不动点公式牛顿迭代法，程序如下：format long;x=sym('x');f=sym('x^3+2*x^2+10*x-20');df=diff(f,x);FX=x-f/df;Fx=inline(FX);disp('x=');x1=0.5;disp(x1);i=0;while 1x0=x1;i=i+1;x1=feval(Fx,x1);disp(x1);if abs(x1-x0)<Epsbreak;endendi,x1运行结果：x=1.0001.6371.3961.4411.137i =4x1 =1.137比较三种方法，牛顿法的收敛性比较好，相比不动点迭代法要构造出收敛的公式比较难，牛顿法迭代次数也较少，收敛速度快，只是对初值的要求很高，几种方法各有利弊，具体采用哪种也需因题而异。