信计091
要求：1、用Matlab语言或你熟悉的其他算法语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。 一、下表给出了飞行中鸭子的上部形状的节点数据，试用三次样条插值函数（自然边界条件）和20次Lagrange插值多项式对数据进行插值。用图示出给定的数据，以及s(x)和L20(x)。
x y x y 0.9 1.3 7.0 2.3 1.3 1.5 8.0 2.25 1.9 1.85 9.2 1.95 2.1 2.1 10.5 1.4 2.6 2.6 11.3 0.9 3.0 2.7 11.6 0.7 3.9 2.4 12 0.6 4.4 2.15 12.6 0.5 4.7 2.05 13.0 0.4 5.0 2.1 13.3 0.25 6.0 2.25 解：>> x=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12 12.6 13.0 13.3];
>> y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25]; （1）三次样条插值法 在MATLAB中编写m文件
function[f,f0]=scyt(x,y,y2_1,y2_N,x0) %y2_1和y2_N分别为自然边界条件
%插值点x的坐标：x0
syms t;
f=0.0;f0=0.0;
if(length(x)==length(y)) n=length(x); else
disp('x和y的维数不相等'); return; end
for i=1:n
if(x(i)<=x0)&&(x(i+1)>=x0) index=i; break; end end
A=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); for i=2:n-1
u(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); end
c(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函数为追赶法，其程序为 function x=zgf(A,b) n = rank(A); for(i=1:n)
if(A(i,i)==0)
disp('Error: 对角有元素为0！'); return; end end;
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); end
d(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); end
x(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)^2 f0 =
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)^2
>> xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline'); >> title('试验一--三次样条插值图示')
试验一--三次样条插值图示32.521.510.5002468101214
（2）用拉格朗日法插值 %定义Lagrange程序
function f=Language(x,y,x0) syms t;
if(length(x)==length(y)) n=length(x); else
disp('x和y的维数不相等');
return; end f=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); else
f=collect(f); f=vpa(f,6); end end end
>> Language(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
?10?7二、已知Wilson矩阵A???8??7解x??1111?。
T77??32??23?565??，且向量b???，则方程组Ax?b有准确?33?6109????5910??31?8⑴用Matlab 内部函数求A，A的所有特征值和cond(A)2；
78.17.2??10?7.085.04?65?，解方程组(A??A)(x??x)?b，并求出向量⑵令A??A???85.989.899???6.994.9999.98???x和?x2，从理论结果和实际计算结果两方面分析方程组Ax?b解的相对误差
?x2x2与A的相对误差?A2A2的关系；
⑶再改变扰动矩阵?A（其元素的绝对值不超过0.005），重复第2问。 解：（1）A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; b=[32;23;33;31]; M=det(A) M = 1
A的所以特征值： >> D=eig(A) D =
0.0102 0.8431 3.8581 30.2887
条件数 >> n=30.2887/0.0102 n =
2.9695e+003
所以A的行列式为1，cond(A)2=2.9695e+003
(2) >> B=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 4.99 9 9.98]; >> b=[32;23;33;31];
>> [rank(B),rank([B,b])] ans = 4 4 >> x=B\\b x =-5.4761 11.4934 -1.4292
2.4838
求?x 假设X= ?x
>> x=B\\b;x1=[1;1;1;1];X=x-x1 X =
-6.4761 10.4934 -2.4292
1.4838 求 ?x2 >>norm(X) ans =
12.6552 12.655就是?x2。 %求?x2x2