数值完成

数值分析上机题目及参考答案一、下表给出了飞行中鸭子的上部形状的节点数据，试用三次样条插值函数（自然边界条件）和20次Lagrange插值多项式对数据进行插值。用图示出给定的数据，以及()sx和20()Lx。x0.91.31.92.12.63.03.94.44.75.06.0y1.31.51.852.12.62.72.42.152.052.12.25x7.08.09.210.511.311.61212.613.013.3y2.32.251.951.40.90.70.60.50.40.25编写三次样条函数m文件：function[f,f0]=scyt(x,y,y2_1,y2_N,x0)%y2_1和y2_N分别为自然边界条件%插值点x的坐标：x0symst;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等');return;endfori=1:nif(x(i)=x0)&&(x(i+1)=x0)index=i;break;endendA=diag(2*ones(1,n));A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);fori=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);endc(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;m=zgf(A,c);h=x(index+1)-x(index);f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h;f0=subs(f,'t',x0);其中的zgf函数为追赶法，其程序为functionx=zgf(A,b)n=rank(A);for(i=1:n)if(A(i,i)==0)disp('Error:对角有元素为0！');return;endend;d=ones(n,1);a=ones(n-1,1);c=ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1)=A(n,n);for(i=2:n)d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1)=b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1)=(b(i,1)-c(i,1)*x(i+1,1))/d(i,1);end在matlab命令窗口输入：x=[0.91.31.92.12.63.03.94.44.75.06.07.08.09.210.511.311.61212.613.013.3];y=[1.31.51.852.12.62.72.42.152.052.12.252.32.251.951.40.90.70.60.50.40.25];[f,f0]=scyt(x,y,0,0,x0,3.5)%其中给出了x=3.5处y的值f=1000/729*(27/5*t-1377/100)*(t-39/10)^2+1000/729*(522/25-24/5*t)*(t-3)^2+100/81*(-6396162352027119/288230376151711744*t+19188487056081357/288230376151711744)*(39/10-t)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*t)*(t-3)^2f0=2.5851三次样条插值图像为plot(x,y,'r*');grid拉格朗日插值函数程序，定义m文件functionf=lglr(x,y,x0)%插值点x的坐标：x0%求得的拉格朗日插值多项式在x0处的值，fsymst;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等');return;endf=0.0;for(i=1:n)l=y(i);for(j=1:i-1);l=l*(t-x(j))/(x(i)-x(j));end;f=f+l;simplify(f);if(i==n)if(nargin==3)f=subs(f,'t',x0)elsef=collect(f);f=vpa(f,6);%插值多项式的系数化为6位小数endendend在命令窗口输入f=lglr(x,y)f=-337.863*t^2+219.386*t-.171512e-3*t^13+.706734*t^8+.152386e-4*t^14-.105677e-5*t^15-.271334*t^9+304.616*t^3-.108815e-1*t^11+.622166e-1*t^10+65.2964*t^5+.152863e-2*t^12-13.7172*t^6+.477508*t^7+.559234e-7*t^16+.769797e-14*t^20-.985682e-12*t^19-60.7814-176.003*t^4+.588899e-10*t^18-.217905e-8*t^17%拉格朗日20次公式f=lglr(x,y,3.5)f=194.2495分析结果知道，拉格朗日插值次数达到了20次，但是和三次样条的结果相差很多，拉格朗日的插值次数不是越高越好。二、已知Wilson矩阵1078775658610975910A，且向量32233331b，则方程组Axb有准确解1111Tx。⑴用Matlab内部函数求A，A的所有特征值和2()condA；A:A=[10787;7565;86109;75910];B=det(A)B=1A的所以特征值：D=eig(A)D=0.01020.84313.858130.2887因为原矩阵是对称的，所以条件数：30.2887/0.0102ans=2.9695e+003条件数Cond（A）=Y1/Yn=30.2887/0.0102=2.9695e+003⑵令1078.17.27.085.046585.989.8996.994.9999.98AA，解方程组()()AAxxb，并求出向量x和2x，从理论结果和实际计算结果两方面分析方程组Axb解的相对误差22xx与A的相对误差22AA的关系；A=[10787;7565;86109;75910];a=[000.10.2;0.080.0400;0-0.02-0.0110;-0.01-0.010-0.02]B=a+A;C=[Bb]b=[32;23;33;31]rank(B)ans=4rank(B)=rank(C)所以扰动矩阵有唯一解rank(C)ans=4x=B\bx=-5.476111.4934-1.42922.4838x=B\b;x1=[1;1;1;1];X=x-x1%求x（设X=x）X=-6.476110.4934-2.42921.4838norm(X)%求2xans=12.655212.655就是2x。norm(X)/norm(x1)%求22xxans=6.3276所以22xx=6.3276norm(a)ans=0.2244norm(A)ans=30.2887norm(a)/norm(A)%求22AAans=0.0074所以22AA=0.0074inv(A)%求A的逆矩阵，下求||)(||1||||||||||||11AAxAAd=inv(A);norm(d)*norm(a)*norm(x)ans=288.48581-norm(d*(a))ans=-12.3693288.4858/-12.3693ans=-23.3227所以||)(||1||||||||||||11AAxAA=-23.3227norm(X)ans=0.0074所以2x||)(||1||||||||||||11AAxAA⑶再改变扰动矩阵A（其元素的绝对值不超过0.005），重复第2问。改变A，取a2=[00.0020.0010.003;0.0010.00400.001;0.003-0.004-0.0010;-0.001-0.0020-0.003]B2=a2+A;C2=[Bb]b=[32;23;33;31]rank(B2)ans=4rank(C2)ans=4rank(B2)=rank(C2)所以扰动矩阵有唯一解x2=B2\bx2=1.06490.88931.02720.9859x=B\b;x1=[1;1;1;1];X=x-x1%求x（设X=x）X2=-6.476110.4934-2.42921.4838norm(X2)%求2xnorm(X2)ans=12.655212.6552就是2xnorm(X)/norm(x1)%求22xxnorm(X2)/norm(x1)ans=6.3276所以22xx=6.3276norm(a2)ans=0.0071norm(a)/norm(A)%求22AAnorm(a2)/norm(A)ans=2.3336e-004所以22AA=2.3336e-004norm(d)*norm(a2)*norm(x2)ans=9.08751-norm(d*(a2))ans=0.6943norm(X2)ans=12.65529.0875/0.6943ans=13.0887所以||)(||1||||||||||||11AAxAA=13.08872x||)(||1||||||||||||11AAxAA三、解三对角线性方程组的追赶法及其应用⑴编写解三对角线性方程组的追赶法的通用程序，并应用于方程组12345210001121000012100001210000120xxxxx，检验程序的正确性；（解为52111,,,,63236Tx）解：由追赶法的递推公式（课本p.160）编写m文件a=[2,2,2,2,2];b=[-1,-1,-1,-1];c=[-1,-1,-1,-1];f=[1,0,0,0,0];n=length(a);b=[0,b];u(1)=f(1)/a(1);v(1)=c(1)/a(1);fork=2:n-1u(k)=(f(k)-u(k-1)*b(k))/(a(k)-v(k-1)*b(k));v(k)=c(k)/(a(k)-v(k-1)*b(k));endu(n)=(r(n)-u(n-1)*b(n))/(a(n)-v(n