用bvp4c Matlab求解边界值

Question

我试图解决一个系统的边界值问题q' = f(q(t), a(t))与输入a在Matlab中使用bvp4c。其中q = [q1, q2, q1_dot, q2_dot]'

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

img: System Equation of the One Pendulum

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: The result should look like this

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)];