搜索
插值法1.下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下symsxl;x1=[01491625364964];y1=[012345678];n=length(x1);Ls=sym(0);fori=1:nl=sym(y1(i));fork=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfork=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/2395008000*x^6(2)三次样条插值，程序如下x1=[01491625364964];y1=[012345678];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=2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x^2+999337332656867/1125899906842624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=√x函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。在[0,1]区间上由上图看不出差别，不妨代入几组数据进行比010203040506070-20020406080100较，取x4=[0:0.2:1]x4=[0:0.2:1];sqrt(x4)%准确值subs(Ls,'x',x4)%拉格朗日插值spline(x1,y1,x4)%三次样条插值运行结果为ans=00.44720.63250.77460.89441.0000ans=00.25040.47300.67060.84551.0000ans=00.24290.46300.66170.84031.0000从这几组数值上可以看出在[0,1]区间上，拉格朗日插值更精确。数据拟合和最佳平方逼近2.有实验给出数据表x0.00.10.20.30.50.81.0y1.00.410.500.610.912.022.46试求3次、4次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。解：(1)三次拟合，程序如下symx;x0=[0.00.10.20.30.50.81.0];y0=[1.00.410.500.610.912.022.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,'--');holdon;plot(x0,y0,'x');xlabel('x');ylabel('y');运行结果S3=-7455778416425083/1125899906842624+1803512222945437/140737488355328*x-655705280524945/140737488355328*x^2+4172976892199509/4503599627370496*x^3图像如下(2)4次多项式拟合symx;x0=[0.00.10.20.30.50.81.0];y0=[1.00.410.500.610.912.022.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');holdon;plot(x0,y0,'x');xlabel('x');ylabel('y');运行结果S3=00.10.20.30.40.50.60.70.80.9100.511.522.5xy3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931990070637/281474976710656*x^2-5965836931688425/1125899906842624*x^3+8491275464650307/9007199254740992*x^4图像如下(3)另一个拟合曲线,新建一个M-file程序如下：function[C,L]=lagran(x,y)w=length(x);n=w-1;L=zeros(w,w);fork=1:n+1V=1;forj=1:n+1ifk~=jV=conv(V,poly(x(j)))/(x(k)-x(j));endendL(k,:)=V;end00.10.20.30.40.50.60.70.80.9100.511.522.5xyC=y*L在命令窗口中输入以下的命令：x=[0.00.10.20.30.50.81.0];y=[1.00.410.500.610.912.022.46];cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x));yy=polyval(cc,xx);plot(xx,yy,'r');holdon;plot(x,y,'x');xlabel('x');ylabel('y');x=[0.00.10.20.30.50.81.0];y=[0.940.580.470.521.002.002.46];%y中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据[C,L]=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);holdon;plot(xx,yy,'b',x,y,'.');图像如下00.10.20.30.40.50.60.70.80.9100.511.522.5xy解线性方程组的直接解法3.线性方程组Ax=b的A及b为A=[1078775658610975910]，b=[32233331]，则解x=[1111].用MATLAB内部函数求detA及A的所有特征值和cond(A)2.若令A+δA=[1078.17.27.085.046585.989.8996.99599.98]，求解(A+δA)(x+δx)=b,输出向量x和||δx||2.从理论结果和实际计算两方面分析线性方程组Ax=b解得相对误差||δx||2/||x||2及A的相对误差||δA||2/||A||2的关系.解：(1)程序如下clear;A=[10787;7565;86109;75910];det(A)cond(A,2)eig(A)输出结果为ans=1ans=2.9841e+003ans=0.01020.84313.858130.2887(2)程序如下A=[1078.17.2;7.085.0465;85.989.899;6.99599.98];b=[32233331]';x0=[1111]';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=[1078.17.2;7.085.0465;85.989.899;6.99599.98];A0=[10787;7565;86109;75910];b=[32233331]';x0=[1111]';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)');gridon;可以估计到方程的根在区间（0,1）;选取迭代初值为x0=0.5；构造不动点迭代公式x(k+1)=(x(k)^2+2-e^x(k))/3;ψ(x)=(x^2+2-e^x)/3;当0x1时，0ψ(x)1;ψ’(x)`1;因此迭代公式收敛。-6-4-20246-200-150-100-50050xx2-3x+2-exp(x)程序如下：formatlong;f=inline('(x^2+2-exp(x))/3')disp('x=');x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while1x0=x;i=i+1;x=feval(f,x);disp(x);ifabs(x-x0)Epsbreak;endendi,x运行结果为f=Inlinefunction:f(x)=(x^2+2-exp(x))/3x=0.200426243099960.272749065098370.253607156584130.258550376264940.257265636335090.257598985162190.257512454514830.257534913615250.257529084167960.257530597238330.257530204510460.257530306445640.257530279987670.25753028685501i=14x=0.25753028685501斯特芬森加速法，程序如下:formatlong;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;while1x0=x;i=i+1;x=feval(f,x);disp(x);ifabs(x-x0)Epsbreak;endendi,x运行结果为x=0.258684427565790.257530317719810.257530285439860.2