How can I solve a 2nd order ODE system with BVP4 library?

Hello!

I'm developing a system and I need to solve a 2nd order ODE system. So, I have a system with 2 equations and 2 variables. First, I have converted the 2nd order variables to 1st order variables. Actually, I'm able to solve the system, but the answer isn't right. That's because the problem has initial conditions and final conditions, and I can't take this in account in the model.

The ODE45 needs only inital conditions. So, I'm using BVP4, and I'm trying to include the right conditions. So, I'll share my code with you:

function bvp4
      options = bvpset('stats','on');

      xlow = 0;
      xhigh = 16;
      solinit = bvpinit(linspace(xlow, xhigh, 50),[0 0 0]);
      sol = bvp4c(@bvp4ode,@bvp4bc,solinit,options);
      xint = linspace(xlow,xhigh);
      Sxint = deval(sol,xint);
      figure
      plot(xint,Sxint(1,:))
      figure
      plot(xint,Sxint(2,:))
      figure
      plot(xint,Sxint(3,:))

% ----------------------------------------------------------

function dx = bvp4ode(t,x)

c_alfa_f = 2149.8183690099601751977414821455/2;
c_alfa_r = 2079.3207560832612588593410087283/2;
M_total = 1.9071e+03;
v_x = (8217*heaviside(t - 231/20))/1000 - (8217*heaviside(t - 1125903425279833/70368744177664))/1000 - (837*heaviside(t - 443/100))/50 + (837*heaviside(t + 4722366482869645/9444732965739290427392))/50 + (heaviside(t - 231/20) - heaviside(t - 443/100))*((775238149725083*t^4)/72057594037927936 - (3253530043737061*t^3)/9007199254740992 + (2374146471846479*t^2)/562949953421312 - (5278954075683331*t)/281474976710656 + 6136806631930043/562949953421312);
d_v_x = (8217*dirac(t - 231/20))/1000 - (8217*dirac(t - 1125903425279833/70368744177664))/1000 - (837*dirac(t - 443/100))/50 + (837*dirac(t + 4722366482869645/9444732965739290427392))/50 + (dirac(t - 231/20) - dirac(t - 443/100))*((775238149725083*t^4)/72057594037927936 - (3253530043737061*t^3)/9007199254740992 + (2374146471846479*t^2)/562949953421312 - (5278954075683331*t)/281474976710656 + 6136806631930043/562949953421312) + (heaviside(t - 231/20) - heaviside(t - 443/100))*((775238149725083*t^3)/18014398509481984 - (9760590131211183*t^2)/9007199254740992 + (2374146471846479*t)/281474976710656 - 5278954075683331/281474976710656);
Iz = 3.8333e+03;
a = 1.2160;
b = 1.5020;
angulo_de_estercamento_total = (652*heaviside(t - 673/120))/37 - (652*heaviside(t - 1033/120))/37 - (652*heaviside(t - 1603/360))/37 + (652*heaviside(t - 1693/360))/37 - ((1440*t)/37 - 5760/37)*(heaviside(t - 4) - heaviside(t - 1603/360)) - ((1440*t)/37 - 7424/37)*(heaviside(t - 673/120) - heaviside(t - 1693/360)) + ((13040*t)/2479 - 467810/7437)*(heaviside(t - 287/24) - heaviside(t - 1033/120))

yaw  = x(1);
beta = x(2);
v    = x(3);

dx = [ v;
       ((-c_alfa_f-c_alfa_r-(M_total*d_v_x))/(M_total*v_x))*beta + ((((-a*c_alfa_f)+(b*c_alfa_r))/(M_total*(v_x.^2)))-1)*v + ((c_alfa_f/(M_total*v_x))*angulo_de_estercamento_total);
       (((-a*c_alfa_f)+(b*c_alfa_r))/Iz)*beta                    + ((-((a^2)*c_alfa_f)-((b^2)*c_alfa_r))/(Iz*v_x))*v       + ((a*c_alfa_f/Iz)*angulo_de_estercamento_total);
       ];

% -------------------------------------------------------------

function res = bvp4bc(x_limite_inferior,x_limite_superior)
res = [ 
        x_limite_inferior(1) % No limite inferior (t = 0s), x(1) = yaw  = 0
        x_limite_inferior(2) % No limite inferior (t = 0s), x(2) = beta = 0
        x_limite_inferior(3) % No limite inferior (t = 0s), x(3) = v    = 0
        ];

So, I have three functions and three variables. The problem is: I have 6 boundary conditions, and I'm not able to include that conditions in BVP4 function.

The boundary conditions are:

          x_limite_inferior(1)+0 % No limite inferior (t = 0s), x(1) = yaw  = 0
          x_limite_inferior(2)+0 % No limite inferior (t = 0s), x(2) = beta = 0
          x_limite_inferior(3)+0 % No limite inferior (t = 0s), x(3) = v    = 0
          x_limite_superior(1)+0 % No limite superior (t = 16s), x(1) = yaw  = 0
          x_limite_superior(2)+0 % No limite superior (t = 16s), x(2) = beta = 0
          x_limite_superior(3)+0 % No limite superior (t = 16s), x(3) = v    = 0

It seems like if I have three variables I'm only able to include 3 boundary conditions. Any one?

Thanks.

[Torsten]

Forget any numerical solver (like bvp4c or ode45) if your problem contains heaviside and dirac functions.

Best wishes

Torsten.

[Lucas Zago]

Thanks for helping Torsten. Is there some method to handle that kind of function? The problem is I have a lot of functions with discontinuities, then the heaviside and dirac function helped me with that. Maybe, there's a way to solve piecewise 2nd order ODE systems functions plus several boundary conditions, but I'm not sure how.