内容发布更新时间 : 2024/6/26 20:11:55星期一 下面是文章的全部内容请认真阅读。

信计08级数值方法计算实习题

要求：1、用Matlab语言或你熟悉的其他算法语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。 一、下表给出了飞行中鸭子的上部形状的节点数据，试用三次样条插值函数（自然边界条件）和20次Lagrange插值多项式对数据进行插值。用图示出给定的数据，以及s(x)和L20(x)。

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 y 编写三次样条函数m文件：function[f,f0]=scyt(x,y,y2_1,y2_N,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=followup(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);

其中的followup函数为追赶法，其程序为 function x=followup(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命令窗口输入：

>> 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];

>> [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)^2 f0 =

2.5851

三次样条插值图像为>> plot(x,y,'r*');grid

拉格朗日插值函数程序，定义m文件

function f=lglr(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; f=f+l; simplify(f); if(i==n)