clear clc
x=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3];
y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25]; n=length(x); for i=1:n-1
h(i)=x(i+1)-x(i); end
for i=1:n-2
k(i)=h(i+1)/(h(i)+h(i+1)); u(i)=h(i)/(h(i)+h(i+1)); end
for i=1:n-2
gl(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1)+k(i)*(y(i+1)-y(i))/h(i)); end
g0=3*(y(2)-y(1))/h(1); g00=3*(y(n)-y(n-1))/h(n-1); g=[g0 gl g00]; g=transpose(g); k1=[k 1]; u1=[1 u];
Q=2*eye(21)+diag(u1,1)+diag(k1,-1); m=transpose(Q\\g); syms X; for i=1:n-1
p1(i)=(1+2*(X-x(i))/h(i))*((X-x(i+1))/h(i))^2*y(i); p2(i)=(1-2*(X-x(i+1))/h(i))*((X-x(i))/h(i))^2*y(i+1); p3(i)=(X-x(i))*((X-x(i+1))/h(i))^2*m(i); p4(i)=(X-x(i+1))*((X-x(i))/h(i))^2*m(i+1); p(i)=p1(i)+p2(i)+p3(i)+p4(i); p(i)=simplify(p(i)); end s1=p(1) s2=p(2) s3=p(3) s4=p(4) for k=1:20
for z=x(k):0.001:x(k+1) q=eval(subs(p(k),'X','z')); plot(z,q,'b') hold on end
end grid on
title('the curve of a bird') xlabel('x') ylabel('p')
结果
s1 =
- (17844995274534675*X^3)/72057594037927936 + (96362974482487425*X^2)/144115188075855872 - (17917369016919769*X)/288230376151711744 + 2867531086713665431/2882303761517117440
s2 =
(460567398537687325*X^3)/486388759756013568 - (1293838210680474065*X^2)/324259173170675712 + (3887388668203714993*X)/648518346341351424 - 31704300895075962577/19455550390240542720 s3 =
- (13314362933111975*X^3)/4503599627370496 + (164459234764219375*X^2)/9007199254740992 - (653533884670203241*X)/18014398509481984 + 4529382045511240869/180143985094819840 s4 =
- (1005732113240237*X^3)/2251799813685248 + (13776694412914149*X^2)/5629499534213120 -
(692306063940821081*X)/225179981368524800 + 2139601553916769761/1125899906842624000