three1.m

计算方法中的三弯矩方法

clc;%三转角
clear;
format long;
x=[8.125 8.4 9.0 9.485 9.6 9.959 10.166 10.2];
fx=[0.0774 0.099 0.280 0.6 0.708 1.200 1.800 2.177];
f1dao=0.01087;
f8dao=100;
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=8.125:0.005:10.2;
for k=1:55
   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=56:175
   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=176:272
   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=273:295
   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=296:366
   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=367:408
   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=409:416
   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('三转角方程计算（1）中三次样条插值函数')
xlabel ('X轴')
ylabel ('Y轴')
grid