三次、五次多项式插值在工程实践中很常见。求解多项式的系数最直接的方法是根据端点处的约束条件，列出线性方程组，再写成矩阵方程AX=B，然后用通用的方法(如高斯消元法、LU分解等)解矩阵方程。
本博文利用matlab符号计算的功能，给出三次、五次多项式插值的系数解析解(不需要解矩阵方程)，并尽可能减少运算量。
一、三次多项式插值
设三次多项式的表达式：
f
(
u
)
=
a
0
+
a
1
(
u
−
u
s
)
+
a
2
(
u
−
u
s
)
2
+
a
3
(
u
−
u
s
)
3
(1)
f(u)=a_0+a_1(u-u_s)+a_2(u-u_s)^2+a_3(u-u_s)^3 \tag 1
f(u)=a0+a1(u−us)+a2(u−us)2+a3(u−us)3(1)
这里为什么不设三次多项式表达式为
f
(
u
)
=
a
0
+
a
1
u
+
a
2
u
2
+
a
3
u
3
f(u)=a_0+a_1u+a_2u^2+a_3u^3
f(u)=a0+a1u+a2u2+a3u3呢？采用式(1)的表达式，可以大大减少求取多项式系数的计算量，读者可以自行尝试加以对比。
插值的端点条件为：
{
f
(
u
s
)
=
p
s
f
′
(
u
s
)
=
v
s
f
(
u
e
)
=
p
e
f
′
(
u
e
)
=
v
e
(2)
\begin{cases} f(u_s)=p_s\\ f'(u_s)=v_s\\ f(u_e)=p_e\\ f'(u_e)=v_e\\ \tag 2 \end{cases}
⎩
⎨
⎧f(us)=psf′(us)=vsf(ue)=pef′(ue)=ve(2)
利用matlab符号计算功能，解得：
{
a
0
=
p
s
a
1
=
v
s
a
2
=
[
3
(
p
e
−
p
s
)
+
(
2
v
s
+
v
e
)
(
u
s
−
u
e
)
]
/
(
u
e
−
u
s
)
2
a
3
=
−
[
2
(
p
e
−
p
s
)
+
(
v
s
+
v
e
)
(
u
s
−
u
e
)
]
/
(
u
e
−
u
s
)
3
(3)
\begin{cases} a_0=p_s\\ a_1=v_s\\ a_2=[3(p_e-p_s)+(2v_s+v_e)(u_s-u_e)]/(u_e-u_s)^2\\ a_3=-[2(p_e-p_s) + (v_s+v_e)(u_s-u_e)]/(u_e-u_s)^3\\ \tag 3 \end{cases}
⎩
⎨
⎧a0=psa1=vsa2=[3(pe−ps)+(2vs+ve)(us−ue)]/(ue−us)2a3=−[2(pe−ps)+(vs+ve)(us−ue)]/(ue−us)3(3)
二、五次多项式插值
设五次多项式的表达式：
f
(
u
)
=
a
0
+
a
1
(
u
−
u
s
)
+
a
2
(
u
−
u
s
)
2
+
a
3
(
u
−
u
s
)
3
+
a
4
(
u
−
u
s
)
4
+
a
5
(
u
−
u
s
)
5
(4)
f(u)=a_0+a_1(u-u_s)+a_2(u-u_s)^2+a_3(u-u_s)^3+a_4(u-u_s)^4+a_5(u-u_s)^5 \tag 4
f(u)=a0+a1(u−us)+a2(u−us)2+a3(u−us)3+a4(u−us)4+a5(u−us)5(4)
插值的端点条件为：
{
f
(
u
s
)
=
p
s
f
′
(
u
s
)
=
v
s
f
′
′
(
u
s
)
=
a
s
f
(
u
e
)
=
p
e
f
′
(
u
e
)
=
v
e
f
′
′
(
u
e
)
=
a
e
(5)
\begin{cases} f(u_s)=p_s\\ f'(u_s)=v_s\\ f''(u_s)=a_s\\ f(u_e)=p_e\\ f'(u_e)=v_e\\ f''(u_e)=a_e\\ \tag 5 \end{cases}
⎩
⎨
⎧f(us)=psf′(us)=vsf′′(us)=asf(ue)=pef′(ue)=vef′′(ue)=ae(5)
利用matlab符号计算功能，解得：
{
a
0
=
p
s
a
1
=
v
s
a
2
=
a
s
/
2
a
3
=
[
(
20
(
p
e
−
p
s
)
+
(
8
v
e
+
12
v
s
)
(
u
s
−
u
e
)
+
(
a
e
−
3
a
s
)
(
u
e
2
+
u
s
2
)
+
(
6
a
s
−
2
a
e
)
u
e
u
s
)
]
/
[
2
(
u
e
−
u
s
)
3
]
a
4
=
[
−
(
30
(
p
e
−
p
s
)
+
(
14
v
e
+
16
v
s
)
(
u
s
−
u
e
)
+
(
2
a
e
−
3
a
s
)
(
u
e
2
+
u
s
2
)
+
(
6
a
s
−
4
a
e
)
u
e
u
s
)
]
/
[
2
(
u
e
−
u
s
)
4
]
a
5
=
[
(
12
(
p
e
−
p
s
)
+
(
6
v
e
+
6
v
s
)
(
u
s
−
u
e
)
+
(
a
e
−
a
s
)
(
u
e
2
+
u
s
2
)
+
(
2
a
s
−
2
a
e
)
u
e
u
s
)
]
/
[
2
(
u
e
−
u
s
)
5
]
(6)
\begin{cases} a_0=p_s\\ a_1=v_s\\ a_2=a_s/2\\ a_3=[(20(p_e - p_s) + (8v_e + 12v_s)(u_s - u_e) + (a_e - 3a_s)(u_e^2 + u_s^2) + (6a_s - 2a_e)u_eu_s)]/[2(u_e-u_s)^3]\\ a_4=[ -(30(p_e - p_s) + (14v_e + 16v_s)(u_s - u_e) + (2a_e - 3a_s)(u_e^2 + u_s^2) + (6a_s - 4a_e)u_eu_s)]/[2(u_e-u_s)^4]\\ a_5=[(12(p_e - p_s) + (6v_e + 6v_s)(u_s - u_e) + (a_e - a_s)(u_e^2 + u_s^2) + (2a_s - 2a_e)u_eu_s)]/[2(u_e-u_s)^5]\\ \tag 6 \end{cases}
⎩
⎨
⎧a0=psa1=vsa2=as/2a3=[(20(pe−ps)+(8ve+12vs)(us−ue)+(ae−3as)(ue2+us2)+(6as−2ae)ueus)]/[2(ue−us)3]a4=[−(30(pe−ps)+(14ve+16vs)(us−ue)+(2ae−3as)(ue2+us2)+(6as−4ae)ueus)]/[2(ue−us)4]a5=[(12(pe−ps)+(6ve+6vs)(us−ue)+(ae−as)(ue2+us2)+(2as−2ae)ueus)]/[2(ue−us)5](6)
三、matlab代码
%{
Function: solve_polyInp_coes
Description: 求解三次、五次插值多项式的系数
Input: 插值多项式结构体
Output: 三次、五次插值多项式的系数a,状态sta(1表示成功，0表示失败)
Author: Marc Pony(marc_pony@163.com)
%}
function [a, sta] = solve_polyInp_coes(polyInp)
sta = 1;
us = polyInp.us;
ue = polyInp.ue;
ps = polyInp.ps;
pe = polyInp.pe;
vs = polyInp.vs;
ve = polyInp.ve;
as = polyInp.as;
ae = polyInp.ae;
if abs(ue - us) < 1.0e-8
sta = 0;
a = [];
return;
end
if polyInp.order == 3
a = zeros(1, 4);
temp = zeros(1, 2);
temp(1) = 1.0 / (ue - us) / (ue - us);
temp(2) = temp(1) / (ue - us);
a(1) = ps;
a(2) = vs;
a(3) = (3*(pe - ps) + (2*vs + ve)*(us - ue)) * temp(1);
a(4) = -(2*(pe - ps) + (ve + vs)*(us - ue)) * temp(2);
elseif polyInp.order == 5
a = zeros(1, 6);
temp = zeros(1, 5);
temp(1) = 0.5 / (ue - us) / (ue - us) / (ue - us);
temp(2) = temp(1) / (ue - us);
temp(3) = temp(2) / (ue - us);
temp(4) = ue * ue + us * us;
temp(5) = ue * us;
a(1) = ps;
a(2) = vs;
a(3) = 0.5 * as;
a(4) = (20*(pe - ps) + (8*ve + 12*vs)*(us - ue) + (ae - 3*as)*temp(4) + (6*as - 2*ae)*temp(5)) * temp(1);
a(5) = -(30*(pe - ps) + (14*ve + 16*vs)*(us - ue) + (2*ae - 3*as)*temp(4) + (6*as - 4*ae)*temp(5)) * temp(2);
a(6) = (12*(pe - ps) + (6*ve + 6*vs)*(us - ue) + (ae - as)*temp(4) + (2*as - 2*ae)*temp(5)) * temp(3);
else
disp('仅支持3次，5次多项式')
sta = 0;
a = [];
return;
end
end
clc
clear
close all
%% 求解三次多项式系数符号解
syms us ue ps pe vs ve as ae real
%f(u) = a0 + a1*u + a2*u^2 + a3*u^3
%f'(u) = a1 + 2*a2*u + 3*a3*u^2
A1 = [1, us, us^2, us^3
0, 1, 2*us, 3*us^2
1, ue, ue^2, ue^3
0, 1, 2*ue, 3*ue^2
];
B1 = [ps; vs; pe; ve];
a1 = simplify(A1 \ B1)
%f(u) = a0 + a1*(u - us) + a2*(u - us)^2 + a3*(u - us)^3
%f'(u) = a1 + 2*a2*(u - us) + 3*a3*(u - us)^2
A2 = [1, 0, 0, 0
0, 1, 0, 0
1, (ue - us), (ue - us)^2, (ue - us)^3
0, 1, 2*(ue - us), 3*(ue - us)^2
];
B2 = [ps; vs; pe; ve];
a2 = simplify(A2 \ B2)
%% 求解五次多项式系数符号解
%f(u) = a0 + a1*u + a2*u^2 + a3*u^3 + a4*u^4 + a5*u^5
%f'(u) = a1 + 2*a2*u + 3*a3*u^2 + 4*a4*u^3 + 5*a5*u^4
%f''(u) = 2*a2 + 6*a3*u + 12*a4*u^2 + 20*a5*u^3
A3 = [1, us, us^2, us^3, us^4, us^5
0, 1, 2*us, 3*us^2, 4*us^3, 5*us^4
0, 0, 2, 6*us, 12*us^2, 20*us^3
1, ue, ue^2, ue^3, ue^4, ue^5
0, 1, 2*ue, 3*ue^2, 4*ue^3, 5*ue^4
0, 0, 2, 6*ue, 12*ue^2, 20*ue^3];
B3 = [ps; vs; as; pe; ve; ae];
a3 = simplify(A3 \ B3)
%f(u) = a0 + a1*(u - us) + a2*(u - us)^2 + a3*(u - us)^3 + a4*(u - us)^4 + a5*(u - us)^5
%f'(u) = a1 + 2*a2*(u - us) + 3*a3*(u - us)^2 + 4*a4*(u - us)^3 + 5*a5*(u - us)^4
%f''(u) = 2*a2 + 6*a3*(u - us) + 12*a4*(u - us)^2 + 20*a5*(u - us)^3
A4 = [1, 0, 0, 0, 0, 0
0, 1, 0, 0, 0, 0
0, 0, 2, 0, 0, 0
1, (ue - us), (ue - us)^2, (ue - us)^3, (ue - us)^4, (ue - us)^5
0, 1, 2*(ue - us), 3*(ue - us)^2, 4*(ue - us)^3, 5*(ue - us)^4
0, 0, 2, 6*(ue - us), 12*(ue - us)^2, 20*(ue - us)^3];
B4 = [ps; vs; as; pe; ve; ae];
a4 = simplify(A4 \ B4)
%% 三次、五次多项式解析解测试
polyInp = struct();
polyInp.order = 3;
polyInp.us = 1;
polyInp.ue = 5;
polyInp.ps = 3;
polyInp.pe = 7;
polyInp.vs = 2;
polyInp.ve = -1;
polyInp.as = 7;
polyInp.ae = 9;
[a, sta] = solve_polyInp_coes(polyInp);
n = 100;
u = linspace(polyInp.us, polyInp.ue, n);
if polyInp.order == 3
pos = a(1) + a(2) * (u - polyInp.us) + a(3) * (u - polyInp.us).^2 + a(4) * (u - polyInp.us).^3;
vel = a(2) + 2.0 * a(3) * (u - polyInp.us) + 3.0 * a(4) * (u - polyInp.us).^2;
figure
subplot(2, 1, 1)
plot(u, pos)
hold on
plot(polyInp.us, polyInp.ps, 'o')
plot(polyInp.ue, polyInp.pe, 'o')
xlabel('u')
ylabel('pos')
title('三次多项式插值')
subplot(2, 1, 2)
plot(u, vel)
hold on
plot(polyInp.us, polyInp.vs, 'o')
plot(polyInp.ue, polyInp.ve, 'o')
xlabel('u')
ylabel('vel')
elseif polyInp.order == 5
pos = a(1) + a(2) * (u - polyInp.us) + a(3) * (u - polyInp.us).^2 + a(4) * (u - polyInp.us).^3 + a(5) * (u - polyInp.us).^4 + a(6) * (u - polyInp.us).^5;
vel = a(2) + 2.0 * a(3) * (u - polyInp.us) + 3.0 * a(4) * (u - polyInp.us).^2 + 4.0 * a(5) * (u - polyInp.us).^3 + 5.0 * a(6) * (u - polyInp.us).^4;
acc = 2.0 * a(3) + 6.0 * a(4) * (u - polyInp.us) + 12.0 * a(5) * (u - polyInp.us).^2 + 20.0 * a(6) * (u - polyInp.us).^3;
figure
subplot(3, 1, 1)
plot(u, pos)
hold on
plot(polyInp.us, polyInp.ps, 'o')
plot(polyInp.ue, polyInp.pe, 'o')
xlabel('u')
ylabel('pos')
title('五次多项式插值')
subplot(3, 1, 2)
plot(u, vel)
hold on
plot(polyInp.us, polyInp.vs, 'o')
plot(polyInp.ue, polyInp.ve, 'o')
xlabel('u')
ylabel('vel')
subplot(3, 1, 3)
plot(u, acc)
hold on
plot(polyInp.us, polyInp.as, 'o')
plot(polyInp.ue, polyInp.ae, 'o')
xlabel('u')
ylabel('acc')
else
end
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
