数值方法计算实习题1

信计091龚立丽200900901004数值方法计算实习题要求：1、用Matlab语言或你熟悉的其他算法语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。一、下表给出了飞行中鸭子的上部形状的节点数据，试用三次样条插值函数（自然边界条件）和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解：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];（1）三次样条插值法在MATLAB中编写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中输入指令[f,f0]=scyt(x,y,0,0)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得三次样条插值函数S(x)=1000/729*(27/5*x-1377/100)*(x-39/10)^2+1000/729*(522/25-24/5*x)*(x-3)^2+100/81*(-6396162352027119/288230376151711744*x+19188487056081357/288230376151711744)*(39/10-x)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*x)*(x-3)^2xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline');title('试验一--三次样条插值图示')0246810121400.511.522.53试验一--三次样条插值图示（2）用拉格朗日法插值%定义Lagrange程序functionf=Language(x,y,x0)symst;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;for(j=i+1:n)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);endendendLanguage(x,y)ans=52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2二、已知Wilson矩阵1078775658610975910A，且向量32233331b，则方程组Axb有准确解1111Tx。⑴用Matlab内部函数求A，A的所有特征值和2()condA；⑵令1078.17.27.085.046585.989.8996.994.9999.98AA，解方程组()()AAxxb，并求出向量x和2x，从理论结果和实际计算结果两方面分析方程组Axb解的相对误差22xx与A的相对误差22AA的关系；⑶再改变扰动矩阵A（其元素的绝对值不超过0.005），重复第2问。解：（1）A=[10787;7565;86109;75910];b=[32;23;33;31];M=det(A)M=1A的所以特征值：D=eig(A)D=0.01020.84313.858130.2887条件数n=30.2887/0.0102n=2.9695e+003所以A的行列式为1，cond(A)2=2.9695e+003(2)B=[1078.17.2;7.085.0465;85.989.899;6.994.9999.98];b=[32;23;33;31];[rank(B),rank([B,b])]ans=44x=B\bx=-5.476111.4934-1.42922.4838求x假设X=xx=B\b;x1=[1;1;1;1];X=x-x1X=-6.476110.4934-2.42921.4838求2xnorm(X)ans=12.655212.655就是2x。%求22xxnorm(X)/norm(x1)ans=6.32766.3276即为22xxnorm(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.322norm(X)ans=0.0074所以2x||)(||1||||||||||||11AAxAA(3)改变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三、解三对角线性方程组的追赶法及其应用⑴编写解三对角线性方程组的追赶法的通用程序，并应用于方程组12345210001121