turnthree.m

计算方法中的三弯矩方法

clc;
clear;
format long;
x=[5.79997266200454   6.00970053330447   6.56195092941116   7.13117188826633   7.28885670047093   7.89060457126068 8.46123974367823   8.75186063074590]
fx=[-5.69051253227686  -5.86969339062953  -6.16597113194670  -6.28264375084247  -6.28759349831078  -6.19354829641297 -5.91565533140666  -5.67311770545967]
f1dao=tan(atan(0.01087)-pi/4);
f8dao=tan(atan(100)-pi/4);
for j=1:7
    h(j)=x(j+1)-x(j)
    fxj(j)=(fx(j+1)-fx(j))/(x(j+1)-x(j))
end
for j=2:7
    lamda(j-1)=h(j)/(h(j-1)+h(j))
    miu(j-1)=h(j-1)/(h(j-1)+h(j))
end
for j=1:6
    fxjj(j)=(fxj(j+1)-fxj(j))/(h(j)+h(j+1))
end
for j=2:7
    g(j)=6*fxjj(j-1)
end
g(1)=6*(fxj(1)-f1dao)/h(1)
g(8)=6*(f8dao-fxj(7))/h(7)
A=[2 1 0 0 0 0 0 0 
   miu(1) 2 lamda(1) 0 0 0 0 0 
   0 miu(2) 2 lamda(2) 0 0 0 0 
   0 0 miu(3) 2 lamda(3) 0 0 0
   0 0 0 miu(4) 2 lamda(4) 0 0 
   0 0 0 0 miu(5) 2 lamda(5) 0 
   0 0 0 0 0 miu(6) 2 lamda(6)
   0 0 0 0 0 0 1 2 ]
M=inv(A)*g'
t=5.79997266200454:0.005:8.75186063074590;
for k=1:42
   sx(k)=M(1)*((x(2)-t(k))^3)/(6*h(1))+M(2)*((t(k)-x(1))^3)/(6*h(1))+(fx(1)-M(1)*(h(1)^2)/6)*(x(2)-t(k))/h(1)+(fx(2)-M(2)*(h(1)^2)/6)*(t(k)-x(1))/h(1);
end
for k=43:153
   sx(k)=M(2)*((x(3)-t(k))^3)/(6*h(2))+M(3)*((t(k)-x(2))^3)/(6*h(2))+(fx(2)-M(2)*(h(2)^2)/6)*(x(3)-t(k))/h(2)+(fx(3)-M(3)*(h(2)^2)/6)*(t(k)-x(2))/h(2);
end
for k=154:267
   sx(k)=M(3)*((x(4)-t(k))^3)/(6*h(3))+M(4)*((t(k)-x(3))^3)/(6*h(3))+(fx(3)-M(3)*(h(3)^2)/6)*(x(4)-t(k))/h(3)+(fx(4)-M(4)*(h(3)^2)/6)*(t(k)-x(3))/h(3);
end
for k=268:298
   sx(k)=M(4)*((x(5)-t(k))^3)/(6*h(4))+M(5)*((t(k)-x(4))^3)/(6*h(4))+(fx(4)-M(4)*(h(4)^2)/6)*(x(5)-t(k))/h(4)+(fx(5)-M(5)*(h(4)^2)/6)*(t(k)-x(4))/h(4);
end
for k=299:419
   sx(k)=M(5)*((x(6)-t(k))^3)/(6*h(5))+M(6)*((t(k)-x(5))^3)/(6*h(5))+(fx(5)-M(5)*(h(5)^2)/6)*(x(6)-t(k))/h(5)+(fx(6)-M(6)*(h(5)^2)/6)*(t(k)-x(5))/h(5);
end
for k=420:533
   sx(k)=M(6)*((x(7)-t(k))^3)/(6*h(6))+M(7)*((t(k)-x(6))^3)/(6*h(6))+(fx(6)-M(6)*(h(6)^2)/6)*(x(7)-t(k))/h(6)+(fx(7)-M(7)*(h(6)^2)/6)*(t(k)-x(6))/h(6);
end
for k=534:591
   sx(k)=M(7)*((x(8)-t(k))^3)/(6*h(7))+M(8)*((t(k)-x(7))^3)/(6*h(7))+(fx(7)-M(7)*(h(7)^2)/6)*(x(8)-t(k))/h(7)+(fx(8)-M(8)*(h(7)^2)/6)*(t(k)-x(7))/h(7);
end
plot(t,sx)
title('三弯矩方程计算（2）中三次样条插值函数')
xlabel ('X轴')
ylabel ('Y轴')
grid