文章目录

• 一、三次多项式插值
• 二、五次多项式插值
• 三、matlab代码

一、三次多项式插值

  设三次多项表达式：
f(u)=a0+a1(u−us)+a2(u−us)2+a3(u−us)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)=a0+a1u+a2u2+a3u3f(u)=a_0+a_1u+a_2u^2+a_3u^3f(u)=a0​+a1​u+a2​u2+a3​u3呢？采用式(1)的表达式，可以大大减少求取多项式系数的计算量，读者可以自行尝试加以对比。

  插值的端点条件为：
{f(us)=psf′(us)=vsf(ue)=pef′(ue)=ve(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符号计算功能，解得：
{a0=psa1=vsa2=[3(pe−ps)+(2vs+ve)(us−ue)]/(ue−us)2a3=−[2(pe−ps)+(vs+ve)(us−ue)]/(ue−us)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)

二、五次多项式插值

  设五次多项表达式：
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)=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(us)=psf′(us)=vsf′′(us)=asf(ue)=pef′(ue)=vef′′(ue)=ae(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符号计算功能，解得：
{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)\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)

三、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