计算方法实验报告习题1(浙大版)

计算方法实验报告

实验名称： 实验1 从函数表出发进行插值 1 引言

某个实际问题中，函数f (x)在区间[a,b]上存在且连续，但难以找到其表达式，只能通过实验和观测得到有限点上的函数表。有些情况虽然可以写出表达式，但结构复杂，使用不方便。所以希望构造简单函数P (x)作为f (x)的近似值。插值法是解决此类问题的一种方法。

设函数y=在插值区间[a,b]上连续，且在n+1个不同的插值节点a≤x 0,x 1,…,x n ≤b 上分别取值y 0,y 1,…,y n 。目的是要在一个性质优良、便于计算的插值函数类Φ中，求一简单函数P (x)，满足插值条件P (x i )=y i (i=0,1,…,n)，而在其他点x≠x i 上，作为f (x)近似值。求插值函数P (x)的方法称为插值法[1]。

2 实验目的和要求

运用Matlab 编写m 文件，定义三种插值函数，要求一次性输入整张函数表，并利用计算机选择在插值计算中所需的节点。分别通过分段线性插值、分段二次插值和全区间上拉格朗日插值计算f (0.15)，f (0.31)，f (0.47)的近似值。

3 算法原理与流程图

（1）原理 1.线性插值

当给定了n+1个点x 0

1

1

11

1)()(------+--=≈i i i i i i i i x x x x y x x x x y x P x f

这种分段低次插值叫分段线性插值。 2.分段二次插值

当给定了n+1个点x 0

专业：电气工程及其自动化 姓名： 李X

∑∏+-=+≠-

=????

?

????

?

????

??

--=≈11

112)()(i i k i k j i j j k j k x x x x y x P x f 这种分段低次插值叫分段二次插值。 3.全区间上拉格朗日插

对节点x i (i=0,1,…,n)中任一点x k (0≤k≤n)，作一n 次多项式l k (x)，使它在该点上的取值为1，在其余点x i (i=0,1,…,k -1,k+1,…,n)上取值为零。对应于每一节点x k (k=0,1,…,n)，都能写出一个满足此条件的多项式，这样写出n+1个多项式l 0(x),l 1(x),…,l n (x)。

)()()()(1100x l y x l y x l y x P n n n ++=

拉格朗日n 次插值多项式（对于全区间上的插值，n 取函数表的长度）

∑=+-+---------=n

k n k k k k k k n k k k

n x x x x x x x x x x x x x x x x y x L 0110110)

())(()()

())(()()(

（2）流程图

j=2,3,…n-2

P ← 0.5*(x j +x j+1)

x

I ←j

Yes

No

j

输入x i ,y i (i=1,2,…,n )及x

I ← n-1

k=1,2,…n

P ← 1

输入x i ,y i (i=1,2,…,n )及x

按公式计算P 2(x)输出x,P 2(x)≈2()(x

P x f j=2,3,…n-2

P ← x j

x

I ←j Yes

No

j

I ← n-1

按公式计算P 1(x)

输出x,P 1(x)1

1)()(-=≈i y x P x f i=1,2,…n i ≠k

p=p*(x-x 0(i))/(x 0(k)-x 0(i))

Yes

No

i

k

输出S

S=0

S=S+y 0(k)*p

分段线性插值法

分段二次插值法

拉格朗日全区间上的插值法

4 程序代码及注释

1.分段线性插值法

%分段线性插值

function y=fdxx(x0,y0,x) %定义函数

p=length(y0);n=length(x0);m=length(x); %计算函数表和x 的长度

if p~=n error('数据输入有误，请重新输入');

%若函数表的x与y长度不一致则输入有误

else

fprintf('分段线性插值\n\n');

for t=1:m %利用循环计算每个x的插值结果

z=x(t);

if zx0(n)

fprintf('x(%d)超出范围;\n',t);

break;

%若x不在函数表范围内，则插值结果将不准确

end

for i=1:n-1

if z

break;

%选取合适的两点使x(i)

end

end

%注：若x不在函数表范围内，则i=n-1

y(t)=y0(i)*(z-x0(i+1))/(x0(i)-x0(i+1))+y0(i+1)*(z-x0(i))/(x0(i+1)-x0(i));

%按照分段线性插值公式求解y

fprintf('y(%d)=%f\nx1=%.3f y1=%.3f,x2=%.3f

y2=%.3f\n\n',t,y(t),x0(i),y0(i),x0(i+1),y0(i+1));

%输出插值结果和所需的节点

end

end

注：若要在Matlab中直接调用插值函数，需要事先将习题1函数表输入

2.分段二次插值法

%分段二次插值

function y=fd2(x0,y0,x) %定义函数分段二次插值

p=length(y0);n=length(x0);m=length(x); %计算函数表和x的长度

if p~=n error('数据输入有误，请重新输入');

%若函数表的x与y长度不一致则输入有误

else

fprintf('分段二次差值\n\n');

for t=1:m %运用循环求解所有点的插值

z=x(t);

if zx0(n)

fprinf('x(%d)超出范围;\n',t);

break;

%如果x不在函数表范围内无法插值

end

i=n-1; %若下列循环i不变，则i=n-1，

for j=1:n-2

if z<(x0(j+1)+x0(j))/2

i=j;

end

end

%选取i使得x(i-1),x(i),x(i+1)是距x最近的三个点

s=0.0;

for k=i-1:i+1

p=1;

for l=i-1:i+1

if l~=k

p=p*(z-x0(l))/(x0(k)-x0(2));

end

end

s=s+y0(k)*p;

end

%根据分段二次插值公式求y

y(t)=s;

fprintf('y(%d)=%f\nx1=%.3f y1=%.3f\nx2=%.3f y2=%.3f\nx3=%.3f y3=%.3f\n\n',t,y(t),x0(i-1),y0(i-1),x0(i),y0(i),x0(i+1),y0(i+1));

%输出结果和所需节点

end

end

3.拉格朗日全区间上插值

%拉格朗日插值

function y=lagrange(x0,y0,x) %定义函数

p=length(y0);n=length(x0);m=length(x); %计算函数表和x的长度

if p~=n error('数据输入有误，请重新输入');

%若函数表的x与y长度不一致则输入有误

else

fprintf('拉格朗日插值\n');

for t=1:m %利用循环计算每个x的插值

s=0.0;

z=x(t);

for k=1:n

p=1;

for i=1:n

if i~=k

p=p*(z-x0(j))/(x0(k)-x0(i));

end

end

s=s+y0(k)*p;

end%根据拉格朗日插值公式求解y

y(t)=s;

fprintf('y(%d)=%f\n',t,y(t)); %输出插值结果

end

end

5算例分析

x0.00.1 0.195 0.3 0.401 0.5

f(x)0.398940.39695 0.39142 0.38138 0.36812 0.35206 从函数表出发计算f(0.15)，f(0.31)，f(0.47)的近似值

>> x0=[0.0 0.1 0.195 0.3 0.401 0.5];

>> y0=[0.39894 0.39695 0.39142 0.38138 0.36812 0.35206];

>> x=[0.15 0.31 0.47];

%首先在Matlab中输入函数表和待插值变量x

%确保Current Directory 中存放分别含有上述代码的三个m文件

（1）分段线性插值

>> y1=fdxx(x0,y0,x) %调用分段线性插值函数得到下列结果

分段线性插值

y(1)=0.394039

x1=0.100 y1=0.397,x2=0.195 y2=0.391

y(2)=0.380067

x1=0.300 y1=0.381,x2=0.401 y2=0.368

y(3)=0.356927

x1=0.401 y1=0.368,x2=0.500 y2=0.352

y1 =

0.39404 0.38107 0.35693

（2）分段二次插值

>> y=fd2(x0,y0,x) %调用分段二次插值函数得到下列结果

分段二次差值

y(1)=0.394554

x1=0.195 y1=0.391

x2=0.300 y2=0.381

x3=0.401 y3=0.368

y(2)=0.380225

x1=0.195 y1=0.391

x2=0.300 y2=0.381

x3=0.401 y3=0.368

y(3)=0.357247

x1=0.300 y1=0.381

x2=0.401 y2=0.368

x3=0.500 y3=0.352

y =

0.39455 0.38122 0.35725

（3）全区间上拉格朗日插值

>> y=lagrange(x0,y0,x) %调用拉格朗日性插值函数得到下列结果

拉格朗日插值

y(1)=0.394473

y(2)=0.380219

y(3)=0.357222

y =

0.39447 0.38122 0.35722

（4）测试示例

>> x=[0.15 0.31 0.51];

>> y=fdxx(x0,y0,x)

分段线性插值

y(1)=0.394039

x1=0.100 y1=0.397,x2=0.195 y2=0.391

y(2)=0.380067

x1=0.300 y1=0.381,x2=0.401 y2=0.368

x(3)超出范围;

>> y1=[ 0.39894 0.39695 0.39142 0.38138 0.36812 ];

>> y=fdxx(x0,y1,x)

??? Error using ==> fdxx

数据输入有误，请重新输入

6讨论与结论

1.对程序运行的时间加以比较

使用tic,toc函数计算下列四种方法计算上述问题所运行的时间

函数interp1(x0,y0,x) fdxx(x0,y0,x) fd2(x0,y0,x) lagrange(x0,y0,x)

运行时间(s) 0.29699 0.062999 0.076434 0.034817 0.25558 0.072315 0.071896 0.038289 0.24297 0.069953 0.085047 0.048222

从三次实验结果可知，interp1函数所用时间比编写的fdxx函数长（interp1函数为matlab求分段线性插值的函数）；而fdxx和另外两个函数时间比较接近

2.程序优化

由分段线性插值和分段二次插值的原理，x取值在函数表范围内时，插值结果有意义，而当x 取值在函数表范围以外，利用分段线性插值公式仍可以进行运算并得到一个值，但其结果不准确；分段二次插值则无法找到三个合适的点以求插值，不予以输出结果；若输入的函数表x 与y的长度不相等，则无法插值。所以加入判断

p=length(y0);n=length(x0);

if p~=n error('数据输入有误，请重新输入');

if zx0(n)

fprintf('x(%d)超出范围;\n',t);

break;

end

从而避免了因输入错误导致的结果错误。由上述（4）测试示例可看到效果。

3.结果比较

（下方的图为放大图）图中data1为原函数表的数据，data2-data4分别为分段线性插值、分段二次插值和拉格朗日插值的结果，其x取步长为0.01的52个数据，可知插值结果与原数据变化规律接近，在给定范围内均可以作其近似函数。

参考文献

[1] 易大义，沈云宝，李有法. 计算方法(第2版)，浙江大学出版社. pp.29-53.

[2] 廖波. 数值计算与最优化原理——MA TLAB实现，北京邮电大学出版社. pp.57-58.