三次、五次多项式插值在工程实践中很常见。求解多项式的系数最直接的方法是根据端点处的约束条件，列出线性方程组，再写成矩阵方程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)文章来源：https://www.toymoban.com/news/detail-464142.html
三、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
文章来源地址https://www.toymoban.com/news/detail-464142.html
到了这里，关于三次、五次多项式插值(附代码)的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
