分段五次多項式插值（MATLAB實現(xiàn)）

分段五次多項式插值（MATLAB實現(xiàn)）?章?錄?、問題描述??給定n+1個點序列(t,p)ti=0,1,...,n。，利?分段五次多項式插值，使得分段多項式經(jīng)過所有點序列。其中，必須單調(diào)遞增，iii?、推導步驟??起點處?階導數(shù)估計：??終點處?階導數(shù)估計：??中間點處?階導數(shù)估計：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??起點處?階導數(shù)估計：??終點處?階導數(shù)估計：??中間點處?階導數(shù)估計：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，?階導數(shù)：p′(t)=a1+sss4s5s2a2(t?t)+3a3(t?t)2+4a4(t?t)3+5a5(t?t)4，?階導數(shù)：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??容易求得系數(shù)：?????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