多篇数值分析实践报告_matlab

数值分析实验报告

拉格朗日插值法

拉格朗日插值法基本原理：

通过平面上不同两点可以确定一条直线，这就是拉格朗日线性插值问题，对于不在同一条直线的三个点得到的插值多项式则为抛物线。

拉格朗日插值的基多项式（即基函数）为：

n

i x x ••

x •x x l i

j j j

i i i ,,2,1,0,)(0 =--=∏

≠=

有了基函数以后就可以直接构造如下多项式：

∑==n

i i i n x l x f p 0

)

()(

该多项式就是拉格朗日插值法所求得的插值多项式。

拉格朗日插值法算法：

1、根据所给点),(i i y x 的坐标依次写出其差值基函数（用for 循环可以轻易解决）

n

i x x ••

x •x x l i

j j j

i i i ,,2,1,0,)(0 =--=∏

≠=

2、将差值基函数与其对应的点的函数值相乘得：

n i x l x f i i ,,2,1,0),()( =

3、将2中各项累加即得差值多项式：

∑==n

i i i n x l x f p 0

)

()(

拉格朗日插值法程序：

function lagrange(A)%A 为一个只有两行的矩阵，第一行为插值点，第二行为插值点对应的函数值 [m,n]=size(A); f=1;

p=0;%两个用到的变量 syms x

for i=1:n

f=(x-A(1,i))*f;

end

for j=1:n

g(j)=f/(x-A(1,j));%求插值基函数的分母

h(j)=subs(g(j),x,A(1,j));%求插值基函数的分子

s(j)=g(j)/h(j)*A(2,j);%插值基函数

end

for k=1:n

s(j)=collect(s(j));%合并同类项

end

for i=1:n

p=p+s(i);

end

fprintf('拉格朗日插值法可得多项式:')

collect(p) % 可用lagrange([1 3 6 8;4 6 9 12])调试

例：有如下表格中有四个插值点及其对应的函数值，用lagrange插值法写出其三次插值多项式：

x 1 3 6 8

i

y 4 6 9 12

i

解：

在matlab命令窗口输入：

lagrange([1 3 6 8;4 6 9 12])

可得运行结果：

拉格朗日插值法可得多项式:

ans =

x^3/70 - x^2/7 + (97*x)/70 + 96/35

牛顿插值法

牛顿插值法基本原理：

拉格朗日插值多项式的理论在许多方面都有应用，是很不错的一种方法，但就插值问题而言，如果增加一个插值点，原先计算插值的多项式就没有用了，即拉格朗日插值法的继承性很差，牛顿插值法就很好地解决了这个问题，具有很好的继承性。

函数)(x f 的差商定义为：

)(][k k x f x f =

k

k k k k k x x x f x f x x f --=

---111][][],[

k

k k k k k k k k x x x x f x x f x x x f --=

------212112]

,[],[],,[

j

k k k j k k j k k j k j k x x x x f x x f x x x f ---+-+----=

]

,,[],,[],,,[111

在差商的基础上可得牛顿插值多项式：

)

())(](,,,[))(](,,[)](,[][)(11010102100100----++--+-+=n n n x x x x x x x x x f x x x x x x x f x x x x f x f x N

牛顿插值法算法：

1、根据差商的定义求出各阶差商。

2、依次写出：1,,1,0),())((110-=---=-n n x x x x x x n n ω。

3、n ω分别与相应的差商相乘，然后累加即可得：

)

())(](,,,[))(](,,[)](,[][)(11010102100100----++--+-+=n n n x x x x x x x x x f x x x x x x x f x x x x f x f x N

牛顿插值法程序：

function newton(A)%A 为两行矩阵，第一行为插值点，第二行为对应插值点的函数值 [m,n]=size(A);

B=zeros(n-1);

k=1;

for i=1:n-1

for j=i+1:n

Z(j)=(A(2,j)-A(2,j-1))/(A(1,j)-A(1,j-k));%求i阶差商

B(i,j)=Z(j);%记住各阶差商

end

for j=i+1:n

A(2,j)=Z(j);

end

k=k+1;

end

f=A(2,1);

g=1;

syms x

for i=1:n-1

g=g*(x-A(1,i));

f=f+B(i,i+1)*g;%求牛顿插值多项式

end

fprintf('牛顿插值法可得多项式:')

f

fprintf('牛顿插值法得到的多项式合并同类项后为:')

collect(f)% newton([1 3 6 8;4 6 9 12])

例：有如下表格中有四个插值点及其对应的函数值，用牛顿插值法写出其三次插值多项式：

x 1 3 6 8

i

y 4 6 9 12

i

解：

在matlab命令窗口输入：

newton([1 3 6 8;4 6 9 12])