Résoudre la valeur limite avec bvp4c Matlab

Question

J'essaie de résoudre un problème de valeur limite d'un système q' = f(q(t), a(t)) avec une entrée a en utilisant bvp4c dans Matlab. Où q = [q1, q2, q1_dot, q2_dot]'

Mon code Matlab ne fonctionne pas correctement. Quelqu'un sait-il comment résoudre ceci?

img: System Equation of the One Pendulum

a est une fonction d'entrée.

l'état initial: q1(0) = pi, q2(0) = 0, q1_dot(0) = 0, q2_dot(0) = 0.

l'état final: q1(te) = 0, q2(te) = 0, q1_dot(te) = 0, q2_dot(te) = 0.

état initial d'entrée

: a(0) = 0, et l'état final: a(te) = 0.

J'ai au total 10 limites. one_pendulum_bc peut prendre seulement 5 et pas plus et moins.

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

je viens de résoudre avec une fonction de modèle.

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