用BVP4C MATLAB求解边界值

Question

我正在尝试使用MATLAB中的BVP4C使用输入a解决系统q' = f(q(t), a(t))的边界值问题。 q = [q1, q2, q1_dot, q2_dot]'

我的MATLAB代码无法正常工作。有人知道如何解决这个问题吗？

img：一个摆的系统方程

a是一个输入函数。

初始状态： q1(0) = pi, q2(0) = 0, q1_dot(0) = 0, q2_dot(0) = 0.

最终状态： q1(te) = 0, q2(te) = 0, q1_dot(te) = 0, q2_dot(te) = 0.

输入初始状态：a(0) = 0,和终端状态：a(te) = 0。

我总共有10个边界。ONE_PENDULUM_BC只能服用5个，而不需要更多。

function one_pendulum
options = [];    % place holder
solinit = bvpinit(linspace(0,1,1000),[pi, 0, 0, 0],0);
sol = bvp4c(@one_pendulum_ode,@one_pendulum_bc,solinit,options);
t = sol.x;
plot(t, sol.y(1,:))
figure(2)
plot(t, sol.y(2,:))
figure(3)
plot(t, sol.y(3,:))
figure(4)
plot(t, sol.y(4,:))
% --------------------------------------------------------------------------
function dydx = one_pendulum_ode(x,y,a)
% Parameter of the One-Pendulum
m1 = 0.3583;    % weight of the pendulum [kg]
J1 = 0.03799;   % moment of inertia [Nms^2]
a1 = 0.43;      % center of gravity [m]
d1 = 0.006588;  % coefficient of friction [Nms]
g  = 9.81;      % gravity [m/s^2]
dydx = [ y(3)
         y(4)
        (a*a1*m1*cos(y(1)) + a1*g*m1*sin(y(1)) - d1*y(3))/(J1 + a1^2 *m1)
         a ];
% --------------------------------------------------------------------------
function res = one_pendulum_bc(ya,yb,a)
res = [ya(1) - pi
       ya(2)
       ya(4)
       yb(1)
       yb(3)];

img：结果应该看起来像

Answer 1

我只是用模型函数解决了它。

function one_pendulum
options = bvpset('stats', 'on', 'NMax', 3000);

T = 2.5;
sol0 = [
-5.3649
119.5379
-187.3080
91.6670];
solinit = bvpinit(linspace(0, T, T/0.002),[pi, 0, 0, 0], sol0);
sol = bvp4c(@one_pendulum_ode,@one_pendulum_bc,solinit,options, T);
t = sol.x;
ax1 = subplot(2,2,1);
ax2 = subplot(2,2,2);
ax3 = subplot(2,2,3);
ax4 = subplot(2,2,4);
plot(ax1, t, sol.y(1,:))
title(ax1, 'phi');
plot(ax2, t, sol.y(2,:))
title(ax2, 'x');
plot(ax3, t, sol.y(3,:))
title(ax3, 'phi_d');
plot(ax4, t, sol.y(4,:))
title(ax4, 'x_d');
length = size(t)
k = sol.parameters
for i=1:length(2)
   k5 = -(k(1) * sin(pi) + k(2) * sin(2*pi)  + k(3) * sin(3*pi) +    k(4) * sin(4*pi))/sin(5*pi);
x = t(i);
acc(i)=k(1) * sin(pi*x/T) + k(2) * sin(2*pi*x/T)  + k(3) * sin(3*pi*x/T) + k(4) * sin(4*pi*x/T) + k5* sin(5*pi*x/T);
  end
figure(2)
plot(t,acc)
% --------------------------------------------------------------------------
function dydx = one_pendulum_ode(x,y,k,T)
% Parameter of the One-Pendulum
m1 = 0.3583;    % weight of the pendulum [kg]
J1 = 0.03799;   % moment of inertia [Nms^2]
a1 = 0.43;      % center of gravity [m]
d1 = 0.006588;  % coefficient of friction [Nms]
g  = 9.81;      % gravity [m/s^2]
% Ansatzfunktion 
k5 = -(k(1) * sin(pi) + k(2) * sin(2*pi)  + k(3) * sin(3*pi) + k(4) *      sin(4*pi))/sin(5*pi);
a =    k(1) * sin(pi*x/T) + k(2) * sin(2*pi*x/T)  + k(3) *    sin(3*pi*x/T) + k(4) * sin(4*pi*x/T) + k5* sin(5*pi*x/T);
dydx = [ y(3)
     y(4)
     (a*a1*m1*cos(y(1)) + a1*g*m1*sin(y(1)) - d1*y(3))/(J1 + a1^2 *m1)
     a ];
  % --------------------------------------------------------------------------
 function res = one_pendulum_bc(ya,yb,k,T)
 res = [ya(1) - pi
   ya(2)
   ya(3)
   ya(4)
   yb(1)
   yb(2)
   yb(3)
   yb(4)];