计算方法实验+编程代码

计算方法实验报告

1 在区间[-1,1]上分别取n=10,20，用两组等距节点对龙格函数21()125f x x

=+作多项式插值，对每个n 值，分别画出插值函数及()f x 的图形。

解：n=10时：

n=20时：

2 在区间[-1,1]上分别取n=10,20，用两组等距节点对龙格函数21()125f x x

=

+作分段线性插值，对每个n 值，分别画出插值函数及()f x 的图形。 解：n=10:

n=20时：

3 对龙格函数21()125f x x =

+在区间[-1,1]上取21, 0,1,2k x k k n n

=-+= ，n 分别取10,20，试分别求3次、5次最小二乘拟合多项式，打印出此曲线拟合函数，分别画出此拟合函数及()f x 的图形。 3次最小二乘拟合

n=10时：

n=20时：

5次最小二乘拟合n=10时：

n=20时：

4 取点

21

cos,0,1,2

2(1)

k

k

x k n

n

π

+

==

+

，n分别取10,20，对龙格函数

2

1

()

125

f x

x

=

+

作多项式插值，对每个n值，分别画出插值函数及()

f x的图形。解：n=10时：

n=20时：

5 比较上面三组近似函数，说说你的体会。你能在此基础上做进一步的探索吗，比如n如果继续增加下去，结果会如何？

附注：编程语言不限，但用matlab等语言编程时，不得直接调用现成的插值与逼近函数，需要你在我们课堂教学的基础上，编程实现上述算法。

解：用不同方法进行插值，得出的插值函数差异较大。其中，最小二乘拟合的曲线精度较低，多项式插值拟合的曲线精度较高。曲线拟合的精度不仅和拟合方法有关，还和采样点位置的选取、个数有关。如1，4题都是多项式插值，但是第4题的拟合度最高。n值越大，拟合的函数会更加接近原函数。

程序：

1．等距节点多项式插值

clear;

clc;

syms x

n=input('input n=');

x1=linspace(-1,1,n+1);

y1=1./(1.+25*x1.^2);

yy=zeros(1,n+1);

fx=0;

for i=1:n+1;

lga=1;

for j=1:n+1

if j~=i

lga=lga*(x-x1(1,j))/(x1(1,i)-x1(1,j));

else

end

end

fx=fx+y1(1,i)*lga;

end

disp(fx);

x=-1:0.01:1;

plot(x,eval(fx),'rh')

hold on

a=linspace(-1,1,100);

y2=1./(1.+25*a.^2);

plot(a,y2,'b')

legend('插值函数','龙格函数')

2.分段线性插值

function y=div(x0,y0,x,n)

for j=1:n

if (x>=x0(j))&&(x<=x0(j+1))

y=(x-x0(j+1))/(x0(j)-x0(j+1))*y0(j)+(x-x0(j))/(x0(j+1)-x0(j ))*y0(j+1);

else

continue

end

end

end

clear;

clc;

n=input('input n=');

x0=linspace(-1,1,n+1);

yy=zeros(2001,1);

for i=1:1:2001

xx(i,1)=0.001*(i-1)-1;

y=div(x0,1./(1+25*x0.^2), xx(i,1),n);

disp(y);

yy(i,1)=y;

end

plot(xx,yy,'r','LineWidth',1.5)

a=linspace(-1,1,100);

hold on

plot(a,1./(1+25*a.^2),'b','LineWidth',1.5)

legend('插值函数','龙格函数')

3.

①三次最小二乘拟合

clear ;

n=input('input n=');

x0=ones(n+1,1);

x1=zeros(n+1,1);

x2=zeros(n+1,1);

x3=zeros(n+1,1);

y1=zeros(n+1,1);

for k=0:1:n

x1(1+k,1)=-1+2*k/n;

y1(1+k,1)=1./(1+25*x1(1+k,1).^2);

x2(k+1,1)=x1(k+1,1)*x1(k+1,1);

x3(k+1,1)=x1(k+1,1)*x1(k+1,1)*x1(k+1,1); end

A=[x0 x1 x2 x3];

B=A'*A;

C=A'*(y1);

D=B\C;

x=linspace(-1,1,51);

y=D(4,1)*x.^3+D(3,1)*x.^2+D(2,1)*x+D(1,1);

a=linspace(-1,1,101);

y0=1./(1+25*a.^2);

plot(x,y,'rh')

hold on

plot(a,y0,'b','LineWidth',2)

legend('插值函数','龙格函数')

②五次最小二乘拟合

clear ;

n=input('input n=');

x0=ones(n+1,1);

x1=zeros(n+1,1);

x2=zeros(n+1,1);

x3=zeros(n+1,1);

x4=zeros(n+1,1);

x5=zeros(n+1,1);

y1=zeros(n+1,1);

for k=0:1:n

x1(1+k,1)=-1+2*k/n;

y1(1+k,1)=1./(1+25*x1(1+k,1).^2);

x2(k+1,1)=x1(k+1,1)*x1(k+1,1);

x3(k+1,1)=x1(k+1,1)*x1(k+1,1)*x1(k+1,1); x4(k+1,1)=(x1(k+1,1)).^4;