??三次、五次多項式插值在工程實踐中很常見。求解多項式的系數(shù)最直接的方法是根據(jù)端點處的約束條件，列出線性方程組，再寫成矩陣方程AX=B，然后用通用的方法(如高斯消元法、LU分解等)解矩陣方程。
??本博文利用matlab符號計算的功能，給出三次、五次多項式插值的系數(shù)解析解(不需要解矩陣方程)，并盡可能減少運算量。
一、三次多項式插值
??設(shè)三次多項式的表達(dá)式：
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)
??這里為什么不設(shè)三次多項式表達(dá)式為
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?+a1?u+a2?u2+a3?u3呢？采用式(1)的表達(dá)式，可以大大減少求取多項式系數(shù)的計算量，讀者可以自行嘗試加以對比。
??插值的端點條件為：
{
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?)=ps?f′(us?)=vs?f(ue?)=pe?f′(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?=ps?a1?=vs?a2?=[3(pe??ps?)+(2vs?+ve?)(us??ue?)]/(ue??us?)2a3?=?[2(pe??ps?)+(vs?+ve?)(us??ue?)]/(ue??us?)3?(3)
二、五次多項式插值
??設(shè)五次多項式的表達(dá)式：
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?)=ps?f′(us?)=vs?f′′(us?)=as?f(ue?)=pe?f′(ue?)=ve?f′′(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?=ps?a1?=vs?a2?=as?/2a3?=[(20(pe??ps?)+(8ve?+12vs?)(us??ue?)+(ae??3as?)(ue2?+us2?)+(6as??2ae?)ue?us?)]/[2(ue??us?)3]a4?=[?(30(pe??ps?)+(14ve?+16vs?)(us??ue?)+(2ae??3as?)(ue2?+us2?)+(6as??4ae?)ue?us?)]/[2(ue??us?)4]a5?=[(12(pe??ps?)+(6ve?+6vs?)(us??ue?)+(ae??as?)(ue2?+us2?)+(2as??2ae?)ue?us?)]/[2(ue??us?)5]?(6)文章來源：http://www.zghlxwxcb.cn/news/detail-464142.html
三、matlab代碼
%{
Function: solve_polyInp_coes
Description: 求解三次、五次插值多項式的系數(shù)
Input: 插值多項式結(jié)構(gòu)體
Output: 三次、五次插值多項式的系數(shù)a,狀態(tài)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
%% 求解三次多項式系數(shù)符號解
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)
%% 求解五次多項式系數(shù)符號解
%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
文章來源地址http://www.zghlxwxcb.cn/news/detail-464142.html
到了這里，關(guān)于三次、五次多項式插值(附代碼)的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
