分段五次多项式插值（MATLAB实现）

分段五次多项式插值（MATLAB实现）⽂章⽬录⼀、问题描述  给定n+1个点序列(t,p)ti=0,1,...,n。，利⽤分段五次多项式插值，使得分段多项式经过所有点序列。其中，必须单调递增，iii⼆、推导步骤  起点处⼀阶导数估计：  终点处⼀阶导数估计：  中间点处⼀阶导数估计：v0=(p1−p0)/(t1−t0)(1)v=(p−pn−1)/(t−tn−1)nnn(2)v=(d+dk+1)/2kk(3)  其中，d=(p−pk−1)/(t−tk−1)。kkk  起点处⼆阶导数估计：  终点处⼆阶导数估计：  中间点处⼆阶导数估计：acc0=(v1−v0)/(t1−t0)(4)(5)(6)acc=(v−vn−1)/(t−tn−1)nnnacc=(e+ek+1)/2kk  其中，e=(v−vk−1)/(t−tk−1)。kkk

  设五次多项式为:p(t)=a0+a1(t−t)+a2(t−t)2+a3ttattatt(−)3+(−)4+(−)5，⼀阶导数：p′(t)=a1+sss4s5s2a2(t−t)+3a3(t−t)2+4a4(t−t)3+5a5(t−t)4，⼆阶导数：p′′(t)=2a2+6a3(t−t)+12a4(t−t)2+20a5(t−sssssst)3。s  对于每⼀段五次多项式，利⽤端点处约束:⎧⎪⎪⎪⎪p(t)=p⎪ss⎪⎪⎪p′(t)=vss⎨p′′(t)=ass⎪⎪(7)⎪p(t)=p⎪⎪⎪⎪⎪p′(t)=veeeep′′(t)=aee  容易求得系数：⎧⎪⎪⎪⎪a=p⎪0s⎪⎪⎪a=vs1⎨a=a/2s⎪⎪2(8)⎪a=[20h−(8v+12v)T−(3a−a)T2]/(2T3)⎪⎪⎪3esse⎪a=[−30h+(14v+16v)T+(3a−2a)T2]/(2T4)esse4a5=[12h−6(v+v)T−(a−a)T2]/(2T5)sese  其中，h=p−p,T=t−t。sese三、MATLAB代码clc;clear;closeall;%{symststepspevsveasaeTreal;a=[1,0,0,0,0,00,1,0,0,0,00,0,2,0,0,01,T,T^2,T^3,T^4,T^50,1,2*T,3*T^2,4*T^3,5*T^40,0,2,6*T,12*T^2,20*T^3]\[ps;vs;as;pe;ve;ae];a=[simplify(a(1)),simplify(a(2)),simplify(a(3)),simplify(a(4)),simplify(a(5)),simplify(a(6))]'%}t=[0,2,4,8,10]';pos=[10,20,0,30,40]';dt=0.001;n=length(t);v=zeros(n,1);acc=zeros(n,1);v(1)=(pos(2)-pos(1))/(t(2)-t(1));v(n)=(pos(n)-pos(n-1))/(t(n)-t(n-1));fork=2:n-1v(k)=0.5*((pos(k)-pos(k-1))/(t(k)-t(k-1))+(pos(k+1)-pos(k))/(t(k+1)-t(k)));v(k)=0.5*((pos(k)-pos(k-1))/(t(k)-t(k-1))+(pos(k+1)-pos(k))/(t(k+1)-t(k)));endacc(1)=(v(2)-v(1))/(t(2)-t(1));acc(n)=(v(n)-v(n-1))/(t(n)-t(n-1));fork=2:n-1acc(k)=0.5*((v(k)-v(k-1))/(t(k)-t(k-1))+(v(k+1)-v(k))/(t(k+1)-t(k)));endtArray=[];posArray=[];velArray=[];accArray=[];tArray=[tArray;t(1)];posArray=[posArray;pos(1)];velArray=[velArray;v(1)];accArray=[accArray;acc(1)];fori=1:n-1ts=t(i);te=t(i+1);ps=pos(i);pe=pos(i+1);vs=v(i);ve=v(i+1);as=acc(i);ae=acc(i+1);h=pe-ps;T=t(i+1)-t(i);a0=ps;a1=vs;a2=0.5*as;a3=(20.0*h-(8.0*ve+12.0*vs)*T-(3.0*as-ae)*T^2)/(2.0*T^3);a4=(-30.0*h+(14.0*ve+16.0*vs)*T+(3.0*as-2.0*ae)*T^2)/(2.0*T^4);a5=(12.0*h-6.0*(ve+vs)*T-(as-ae)*T^2)/(2.0*T^5);tt=(t(i)+dt:dt:t(i+1))';ifabs(tt(end)-t(i+1))>1.0e-8tt=[tt;t(i+1)];endtArray=[tArray;tt];posArray=[posArray;a0+(tt-ts).*(a1+(tt-ts).*(a2+(tt-ts).*(a3+(tt-ts).*(a4+a5.*(tt-ts)))))];velArray=[velArray;a1+(tt-ts).*(2.0*a2+(tt-ts).*(3.0*a3+(tt-ts).*(4.0*a4+5.0*a5.*(tt-ts))))];accArray=[accArray;2.0*a2+(tt-ts).*(6.0*a3+(tt-ts).*(12.0*a4+20.0*a5.*(tt-ts)))];endfigure(1)subplot(3,1,1)plot(t,pos