三次、五次多项式插值(附代码)
文章目录
- 一、三次多项式插值
- 二、五次多项式插值
- 三、matlab代码
三次、五次多项式插值在工程实践中很常见。求解多项式的系数最直接的方法是根据端点处的约束条件，列出线性方程组，再写成矩阵方程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-8sta = 0;a = [];return;
endif polyInp.order == 3a = 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 == 5a = 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);
elsedisp('仅支持3次，5次多项式')sta = 0;a = [];return;
endend
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^30, 1, 2*us, 3*us^21, ue, ue^2, ue^30, 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, 00, 1, 0, 01, (ue - us), (ue - us)^2, (ue - us)^30, 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^50, 1, 2*us, 3*us^2, 4*us^3, 5*us^40, 0, 2, 6*us, 12*us^2, 20*us^31, ue, ue^2, ue^3, ue^4, ue^50, 1, 2*ue, 3*ue^2, 4*ue^3, 5*ue^40, 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, 00, 1, 0, 0, 0, 00, 0, 2, 0, 0, 01, (ue - us), (ue - us)^2, (ue - us)^3, (ue - us)^4, (ue - us)^50, 1, 2*(ue - us), 3*(ue - us)^2, 4*(ue - us)^3, 5*(ue - us)^40, 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 == 3pos = 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;figuresubplot(2, 1, 1)plot(u, pos)hold onplot(polyInp.us, polyInp.ps, 'o')plot(polyInp.ue, polyInp.pe, 'o')xlabel('u')ylabel('pos')title('三次多项式插值')subplot(2, 1, 2)plot(u, vel)hold onplot(polyInp.us, polyInp.vs, 'o')plot(polyInp.ue, polyInp.ve, 'o')xlabel('u')ylabel('vel')
elseif polyInp.order == 5pos = 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;figuresubplot(3, 1, 1)plot(u, pos)hold onplot(polyInp.us, polyInp.ps, 'o')plot(polyInp.ue, polyInp.pe, 'o')xlabel('u')ylabel('pos')title('五次多项式插值')subplot(3, 1, 2)plot(u, vel)hold onplot(polyInp.us, polyInp.vs, 'o')plot(polyInp.ue, polyInp.ve, 'o')xlabel('u')ylabel('vel')subplot(3, 1, 3)plot(u, acc)hold onplot(polyInp.us, polyInp.as, 'o')plot(polyInp.ue, polyInp.ae, 'o')xlabel('u')ylabel('acc')
elseend
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！