数值分析第二章实习题
第一题
拉格朗日插值
test.m
程序：
function yy=test(x,y,xx)
n=length(x);
m=length(y);
if n~=m
error(‘x和y的维数必须相同’);
return;
end
p=zeros(1,n);
for i=1:n
t=ones(1,n);
for j=[1:i-1,i+1:n]
t(j)=(xx-x(j))/(x(i)-x(j));
end
p(i)=prod(t); %%li(xx),i=1:n,prod(t)表示向量t各元素分量相乘
end
yy=sum(y.*p);
end
运行结果：
>> yy=test(x,y,[0.2])
yy =
0.9800
nalagr.m
程序：
function yy=nalagr(x,y,xx)
m=length(x);n=length(y);
if m~=n,error(‘向量x与y的长度必须保持一致’);end
s=0;
for i=1:n
t=ones(1,length(xx));
for j=[1:i-1,i+1:n]
t=t.*(xx-x(j))/(x(i)-x(j));
end
s=s+t*y(i);
end
yy=s;
运行结果：
>> xx
xx =
0.2000 0.2800 1.0800 1.0000
>> yy=nalagr(x,y,xx)
yy =
0.9800 0.9622 0.2403 0.3800
牛顿插值
程序：
function [F,yy]= newton(X,Y,xx)
% Newton插值函数
% X为已知数据点的x坐标
% Y为已知数据点的y坐标
% xx为插值点的x坐标
% 函数返回F差商表
% yy为各插值点函数值
n=length(X); m=length(xx); % xx=[0.2 0.28 1.0 1.08]
F=zeros(n,n);
F(:,1)=Y’; %矩阵的第一列用于存放插值向量X的函数值向量Y
for j=2:n
for i=j:n
F(i,j)=(F(i,j-1)- F(i-1,j-1))/(X(i)-X(i-j+1)); %求矩阵下三角部分元素值
end
end
% disp(‘A:’)
% A
for t=1:m % m=4
x=xx(t);
s=0.0;
for k=1:n
w=1.0;
for j=1:k-1
w=w*(x-X(j)); %Wn(x)=(x-x1)*(x-x2)……(x-xn) %n=1,2,3,4,5
end
s=s+F(k,k)*w;
end
ss(t)=s;
end
yy=ss;
F=[X’,F]; %将X向量拼进矩阵F的第一列中
end
运行结果：
>> X=0.2:0.2:1;
>> Y=[0.98 0.92 0.81 0.64 0.38];
>> xx=[0.2 0.28 1.0 1.08];
>> [F,yy]=newton(X,Y,xx)
F =
0.2000 0.9800 0 0 0 0
0.4000 0.9200 -0.3000 0 0 0
0.6000 0.8100 -0.5500 -0.6250 0 0
0.8000 0.6400 -0.8500 -0.7500 -0.2083 0
1.0000 0.3800 -1.3000 -1.1250 -0.6250 -0.5208
yy =
0.9800 0.9622 0.3800 0.2403
牛顿四次插值公式为：
$$
P_4\left( x \right) =0.98−0.3\left( x−0.2 \right) −0.625\left( x−0.2 \right) \left( x−0.4 \right) −
$$
$$
0.2083\left( x−0.2 \right) \left( x−0.4 \right) \left( x−0.6 \right) −0.5208\left( x−0.2 \right) \left( x−0.4 \right) \left( x−0.6 \right) \left( x−0.8 \right)
$$
第二题
n=10
程序：
X=-1:0.2:1;
Y=1./(1+25.*X.^2);
n=length(X);
F=zeros(n,n);
F(:,1)=Y’; %矩阵的第一列用于存放插值向量X的函数值向量Y
for j=2:n
for i=j:n
F(i,j)=(F(i,j-1)- F(i-1,j-1))/(X(i)-X(i-j+1)); %求矩阵下三角部分元素值
end
end
F=[X’,F]; %将X向量拼进矩阵F的第一列中
F
$$
P_{10}\left( x \right) =
$$
$$
0.0385+
$$
$$
0.1018\left( x+1 \right) +
$$
$$
0.2602\left( x+1 \right) \left( x+0.8 \right) +
$$
$$
0.7919\left( x+1 \right) \left( x+0.8 \right) \left( x+0.6 \right) +
$$
$$
2.6867\left( x+1 \right) \left( x+0.8 \right) \left( x+0.6 \right) \left( x+0.4 \right) –
$$
$$
6.3631\left( x+1 \right) \left( x+0.8 \right) \left( x+0.6 \right) \left( x+0.4 \right) \left( x+0.2 \right) –
$$
$$
17.6753\left( x+1 \right) \left( x+0.8 \right) \left( x+0.6 \right) \left( x+0.4 \right) \left( x+0.2 \right) x+
$$
$$
84.8416\left( x+1 \right) \left( x+0.8 \right) \left( x+0.6 \right) \left( x+0.4 \right) \left( x+0.2 \right) x\left( x-2 \right) –
$$
$$
167.9157\left( x+1 \right) \left( x+0.8 \right) \left( x+0.6 \right) \left( x+0.4 \right) \left( x+0.2 \right) x\left( x-0.2 \right) \left( x-0.4 \right) +
$$
$$
220.9417\left( x+1 \right) \left( x+0.8 \right) \left( x+0.6 \right) \left( x+0.4 \right) \left( x+0.2 \right) x\left( x-0.2 \right) \left( x-0.4 \right) \left( x-0.6 \right) –
$$
$$
220.9417\left( x+1 \right) \left( x+0.8 \right) \left( x+0.6 \right) \left( x+0.4 \right) \left( x+0.2 \right) x\left( x-0.2 \right) \left( x-0.4 \right) \left( x-0.6 \right) \left( x-0.8 \right)
$$
$$
$$
运行结果：
>> y2=
0.0385+
0.1018.*(x+1)+
0.2602.*(x+1).*(x+0.8)+
0.7919.*(x+1).*(x+0.8).*(x+0.6)+
2.6867.*(x+1).*(x+0.8).*(x+0.6).*(x+0.4)-
6.3631.*(x+1).*(x+0.8).*(x+0.6).*(x+0.4).*(x+0.2)-17.6753.*(x+1).*(x+0.8).*(x+0.6).*(x+0.4).*(x+0.2).*x+
84.8416.*(x+1).*(x+0.8).*(x+0.6).*(x+0.4).*(x+0.2).*x.*(x-0.2)-167.9157.*(x+1).*(x+0.8).*(x+0.6).*(x+0.4).*(x+0.2).*x.*(x-0.2).*(x-0.4)+
220.9417.*(x+1).*(x+0.8).*(x+0.6).*(x+0.4).*(x+0.2).*x.*(x-0.2).*(x-0.4).*(x-0.6)-220.9417.*(x+1).*(x+0.8).*(x+0.6).*(x+0.4).*(x+0.2).*x.*(x-0.2).*(x-0.4).*(x-0.6).*(x-0.8);
>> plot(x,y1,x,y2)
橙色的曲线为10次牛顿插值多项式，可知其在-1，+1处跟原函数图像有较大差异，误差很大，故并非插值次数n越大越好。
n=20
修改程序中的一行代码为X=-1:0.1:1;
插值公式的求法类似，不再赘述，当n=20时，端点处值的变化更加剧烈，误差更大。
第三题
>> x2
x2 =
1 至 17 列
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
18 至 34 列
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
35 至 51 列
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
52 至 65 列
51 52 53 54 55 56 57 58 59 60 61 62 63 64
>> y2
y2 =
1 至 10 列
0 1.0000 1.4142 1.7321 2.0000 2.2361 2.4495 2.6458 2.8284 3.0000
11 至 20 列
3.1623 3.3166 3.4641 3.6056 3.7417 3.8730 4.0000 4.1231 4.2426 4.3589
21 至 30 列
4.4721 4.5826 4.6904 4.7958 4.8990 5.0000 5.0990 5.1962 5.2915 5.3852
31 至 40 列
5.4772 5.5678 5.6569 5.7446 5.8310 5.9161 6.0000 6.0828 6.1644 6.2450
41 至 50 列
6.3246 6.4031 6.4807 6.5574 6.6332 6.7082 6.7823 6.8557 6.9282 7.0000
51 至 60 列
7.0711 7.1414 7.2111 7.2801 7.3485 7.4162 7.4833 7.5498 7.6158 7.6811
61 至 65 列
7.7460 7.8102 7.8740 7.9373 8.0000
>> plot(x2,y2,’r’,x2,y22,’g’)
>> f=@(x2)x2-0.1667.*x2.*(x2-1)+0.1667.*x2.*(x2-1).*(x2-4)-9.9206*10^(-4).*x2.*(x2-1).*(x2-4).*(x2-9)+3.858.*x2.*(x2-1).*(x2-4).*(x2-9).*(x2-16)-1.0522*10^(-6).*x2.*(x2-1).*(x2-4).*(x2-9).*(x2-16).*(x2-25)+2.1198*10^(-8).*x2.*(x2-1).*(x2-4).*(x2-9).*(x2-16).*(x2-25).*(x2-36)-3.2806*10^(-10).*x2.*(x2-1).*(x2-4).*(x2-9).*(x2-16).*(x2-25).*(x2-36).*(x2-49)
f =
@(x2)x2-0.1667.*x2.*(x2-1)+0.1667.*x2.*(x2-1).*(x2-4)-9.9206*10^(-4).*x2.*(x2-1).*(x2-4).*(x2-9)+3.858.*x2.*(x2-1).*(x2-4).*(x2-9).*(x2-16)-1.0522*10^(-6).*x2.*(x2-1).*(x2-4).*(x2-9).*(x2-16).*(x2-25)+2.1198*10^(-8).*x2.*(x2-1).*(x2-4).*(x2-9).*(x2-16).*(x2-25).*(x2-36)-3.2806*10^(-10).*x2.*(x2-1).*(x2-4).*(x2-9).*(x2-16).*(x2-25).*(x2-36).*(x2-49)
>> f(4)
ans =
1.9996
>> f(64)
ans =
2.4640e+09
当取x=2时，y=1.9996，结果误差小，但是当n=64时，值变得非常大，误差巨大。
原文地址：http://www.cnblogs.com/zhangfurong/p/16886116.html
