turnthree1.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=2:7
    g(j-1)=3*(lamda(j-1)*fxj(j-1)+miu(j-1)*fxj(j))
end
A=[2 miu(1) 0 0 0 0 
   lamda(2) 2 miu(2) 0 0 0
   0 lamda(3) 2 miu(3) 0 0
   0 0 lamda(4) 2 miu(4) 0
   0 0 0 lamda(5) 2 miu(5)
   0 0 0 0 lamda(6) 2 ]
g(1)=g(1)-(lamda(1)*f1dao)
g(6)=g(6)-(miu(6)*f8dao)
C=inv(A)*g'
M=[f1dao C(1) C(2) C(3) C(4) C(5) C(6) f8dao]
t=5.79997266200454:0.005:8.75186063074590;
for k=1:42
   sx(k)=((t(k)-x(2))^2)*(h(1)+2*(t(k)-x(1)))*fx(1)/h(1)^3 + (t(k)-x(1))^2*(h(1)+2*(x(2)-t(k)))*fx(2)/h(1)^3  + (t(k)-x(2))^2*(t(k)-x(1))*M(1)/h(1)^2 + (t(k)-x(2))*(t(k)-x(1))^2*M(2)/h(1)^2;
end
for k=43:153
   sx(k)=((t(k)-x(3))^2)*(h(2)+2*(t(k)-x(2)))*fx(2)/h(2)^3 + (t(k)-x(2))^2*(h(2)+2*(x(3)-t(k)))*fx(3)/h(2)^3  + (t(k)-x(3))^2*(t(k)-x(2))*M(2)/h(2)^2 + (t(k)-x(3))*(t(k)-x(2))^2*M(3)/h(2)^2;
end
for k=154:267
   sx(k)=(t(k)-x(4))^2*(h(3)+2*(t(k)-x(3)))*fx(3)/h(3)^3 + (t(k)-x(3))^2*(h(3)+2*(x(4)-t(k)))*fx(4)/h(3)^3  + (t(k)-x(4))^2*(t(k)-x(3))*M(3)/h(3)^2 + (t(k)-x(4))*(t(k)-x(3))^2*M(4)/h(3)^2;
end
for k=268:298
   sx(k)=(t(k)-x(5))^2*(h(4)+2*(t(k)-x(4)))*fx(4)/h(4)^3 + (t(k)-x(4))^2*(h(4)+2*(x(5)-t(k)))*fx(5)/h(4)^3  + (t(k)-x(5))^2*(t(k)-x(4))*M(4)/h(4)^2 + (t(k)-x(5))*(t(k)-x(4))^2*M(5)/h(4)^2;
end
for k=299:419
   sx(k)=(t(k)-x(6))^2*(h(5)+2*(t(k)-x(5)))*fx(5)/h(5)^3 + (t(k)-x(5))^2*(h(5)+2*(x(6)-t(k)))*fx(6)/h(5)^3  + (t(k)-x(6))^2*(t(k)-x(5))*M(5)/h(5)^2 + (t(k)-x(6))*(t(k)-x(5))^2*M(6)/h(5)^2;
end
for k=420:533
   sx(k)=(t(k)-x(7))^2*(h(6)+2*(t(k)-x(6)))*fx(6)/h(6)^3 + (t(k)-x(6))^2*(h(6)+2*(x(7)-t(k)))*fx(7)/h(6)^3  + (t(k)-x(7))^2*(t(k)-x(6))*M(6)/h(6)^2 + (t(k)-x(7))*(t(k)-x(6))^2*M(7)/h(6)^2;
end
for k=534:591
   sx(k)=((t(k)-x(8))^2)*(h(7)+2*(t(k)-x(7)))*fx(7)/(h(7)^3) + ((t(k)-x(7))^2)*(h(7)+2*(x(8)-t(k)))*fx(8)/h(7)^3  + (t(k)-x(8))^2*(t(k)-x(7))*M(7)/h(7)^2 + (t(k)-x(8))*(t(k)-x(7))^2*M(8)/h(7)^2;
end
plot(t,sx)
title('三转角方程计算（2）中三次样条插值函数')
xlabel ('X轴')
ylabel ('Y轴')
grid