How do I solve this 2nd order ODE with time-dependent parameters?

 2nd ODE to solve
   dx/dt*(a + f(t))* d2x/dt2 + 0.5*b*(dx/dt)^2 + k(t)*(dx/dt)^2 - g(t) = 0
   Boundary condition: dx/dt(0) = v0
 where
   - 't' is the time, 
   - 'x' is the position
   - 'dx/dt' is the velocity
   - 'd2x/dt2' is the acceleration
   - 'a', 'b', 'v0' are constants
   -  'f(t)', 'k(t)' and 'h(t)' are KNOWN functions dependent on 't' 
             (I do now write them because they are quite big)

Since I am new in Matlab, I would like to know how to code this with ode45. Thank you in advance

[Subrata Chakrabarti]

I saw the reply from Torsten which was very nice indeed. I only have a small clarification to ask. What if the force term f(t), which is also a time dependent variable, is obtained from a separate algebraic expression which is running on a for loop with a counter similar to the time step.

Answer 1

[Torsten]

v0 = ...;
tstart = ...;
tend = ...;
a= ...;
b =...;
fun=@(t,y)(g(t)-(0.5*b+k(t))*y(1)^2)/(y(1)*(a+f(t)));
[T,Y]=ode45(fun,[tstart tend],v0);

Best wishes

Torsten.

[Jaime R]

Thank you very much Torsten.

Now I have an error regarding the inputs of the ode45. The time-dependent parameters 'f(t)', 'k(t)' and 'h(t)' are defined as a result of some previous differentations with the time 't' considered as a symbolic variable, so this can't be changed. So, for example, in the end g(t) is defined as:

 g(t) = a + b*t;

but 'a' and 'b' here are huge numbers, like a = 1e+04 and b = 1e+14. So I have the following error:

 Error using odearguments (line 110)
Inputs must be floats, namely single or double.

 Error in ode45 (line 113)
 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

 Error in myFileXX (line XX)
                  [T,Y = ode45(fun, [0 1], v0);

What can I do to solve that? I read something about using:

 g(t) = vpa(g(t));
 f(t) = vpa(f(t));
 k(t) = vpa(k(t));

but it didn't work, the same error still shows up, although these parameters are simplified by that.

Thank you in advance again.

[Torsten]

"matlabFunction" converts symbolic expressions into function handles.

If this doesn't help, we'll have to see how you exactly get g(t),f(t) and k(t).

Best wishes

Torsten.

[Jaime R]

It seems that with matlabFunction it gives another problem regarding the variable y. I show you the full code so the error is clearer.

 syms t        % Variable time
 syms y(t)     % Variable velocity

    %%--- Initial conditions ---
       phi = 12.5e-3;
       v0 = 300;            
       e  = 3e-3;           
       ro = 1580;          
       E  = 43e9;   
       e_r = 0.01466;  
       B = 0.28e-3; 

    %%--- Intermediate calculations ---
       v_T = sqrt(((1 + e_r) * 620e6) /E) - sqrt(E/ro) * e_r;
       R_T = v_T * t;                   
       m_acc = pi * e * ro *(R_T^2);  

       v_L = sqrt (E/ro);           
       R_L = v_L * t;    
       z = 2 * R_L;                   
       E_4_dt = B * ((e_r^2)* B * (0.9^(z/B)-1)) /(log(0.9));

       E_1 = E * e * pi * e_r^2 * (-phi* (phi - 2*v_T*t)) /16; 
       E_2 = pi * R_T^2 * 10e9;        
       E_3 = pi * R_T^2 * 1e6 * e;    

        %%Resolution of the problem  
            g_t = -diff(E_1 + E_2 + E_3, t); 

            % Nonlinear ODE to solve
                  fun = @(t,y)(g_t - (0.5*pi*v_T*e*ro + E_4_dt) * y^2) /(y* (8.33e-3 + m_acc));
                  [T,Y] = ode45(fun, [0 1], v0);

In the end, I would like to plot the results as:

plot(T,Y)

I hope this clarify my doubts. Thank you again for your suggestions Torsten!

[Subrata Chakrabarti]

I saw the reply from Torsten which was very nice indeed. I only have a small clarification to ask. What if any time dependent term within the parantheses like g_t is obtained from a separate algebraic expression which is running on a for loop with a counter similar to the time step.

Answer 2

[Torsten]

Try

syms t y
%%--- Initial conditions ---
       phi = 12.5e-3;
       v0 = 300;            
       e  = 3e-3;           
       ro = 1580;          
       E  = 43e9;   
       e_r = 0.01466;  
       B = 0.28e-3; 
    %%--- Intermediate calculations ---
       v_T = sqrt(((1 + e_r) * 620e6) /E) - sqrt(E/ro) * e_r;
       R_T = v_T * t;                   
       m_acc = pi * e * ro *(R_T^2);  
       v_L = sqrt (E/ro);           
       R_L = v_L * t;    
       z = 2 * R_L;                   
       E_TF_dt = B * ((e_r^2)* B * (0.9^(z/B)-1)) /(log(0.9));
       E_ED = E * e * pi * e_r^2 * (-phi* (phi - 2*v_T*t)) /16; 
       E_DL = pi * R_T^2 * 10e9;        
       E_MC = pi * R_T^2 * 1e6 * e;    
        %%Resolution of the problem  
            g_t = -diff(E_ED + E_DL + E_MC, t); 
f(t,y)=(g_t - (0.5*pi*v_T*e*ro + E_TF_dt) * y^2) /(y* (8.33e-3 + m_acc));
fun=matlabFunction(f);
[T,Y]=ode45(fun,[0 1], v0);

Best wishes

Torsten.

[Jaime R]

It worked Torsten! Thank you very much :)

I was wondering also, since it's a 2nd order ODE, how could I get the 'x' as a function of the time 't' (numerically)? If:

 y = dx/dt

The definition of 'f(t,y)' should change a bit in order to the get 'y' and 'x' with respect to 't', shouldn't it? How could this be coded in Matlab then?

Thank you a lot again Torsten!

[Torsten]

This might work, but I'm not sure:

syms t y(2)
%%--- Initial conditions ---
       phi = 12.5e-3;
       v0 = 300;   
       x0 = ...;         
       e  = 3e-3;           
       ro = 1580;          
       E  = 43e9;   
       e_r = 0.01466;  
       B = 0.28e-3; 
    %%--- Intermediate calculations ---
       v_T = sqrt(((1 + e_r) * 620e6) /E) - sqrt(E/ro) * e_r;
       R_T = v_T * t;                   
       m_acc = pi * e * ro *(R_T^2);  
       v_L = sqrt (E/ro);           
       R_L = v_L * t;    
       z = 2 * R_L;                   
       E_TF_dt = B * ((e_r^2)* B * (0.9^(z/B)-1)) /(log(0.9));
       E_ED = E * e * pi * e_r^2 * (-phi* (phi - 2*v_T*t)) /16; 
       E_DL = pi * R_T^2 * 10e9;        
       E_MC = pi * R_T^2 * 1e6 * e;    
        %%Resolution of the problem  
            g_t = -diff(E_ED + E_DL + E_MC, t); 
f(t,y)=[y(2) ; g_t - (0.5*pi*v_T*e*ro + E_TF_dt) * y(2)^2) /(y(2)* (8.33e-3 + m_acc))];
fun = matlabFunction(f);
[T,Y]=ode45(fun,[0 1], [x0 ; v0]);

Best wishes

Torsten.