ODE45 with very large numbers as constraints

Question

2^nd ODE to solve in MATLAB:

  ( (a + f(t))·d²x/dt² + (b/2 + k(t))·dx/dt ) · dx/dt - 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 not write them because they are quite big)

As an example, using symbolic variables:

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_4 = 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); 

f(t,y)=(g_t - (pi*v_T*e*ro/2 + E_4) * y^2 /(y* (8.33e-3 + m_acc))];
fun=matlabFunction(f);
[T,Y]=ode45(fun,[0 1], v0]);

How can I rewrite this to get x as y=dx/dt? I'm new to Matlab and any help is very welcome !

Answer 1

First, you chould use subs to evaluate a symbolic function. Another approach is to use matlabFunction to convert all symbolic expressions to anonymous functions, as suggested by Horchler.

Second, you're integrating the ODE as if it is 1^st order in dx/dt. If you're interested in x(t) as well as dx/dt(t), then you'll have to modify the function like so:

fun = @(t,y) [y(2);   
             ( subs(g) - (b/2 + subs(k))*y(2)*y(2) ) / ( y(2) * (a + subs(f))) ];

and of course, provide an initial value for x0 = x(0) as well as v0 = dx/dt(0).

Third, the absolute value of the parameters is hardly ever a real concern. IEEE754 double-precision floating point format can effortlessly represent numbers between 2.225073858507201e-308 and 1.797693134862316e+308 (realmin and realmax, respectively). So for the coefficients you gave (O(10¹⁴)), this is absolutely not a problem. You might lose a few digits of precision if you don't take precautions (rescale to [-1 +1], reformulate the problem in different units, ...), but the relative error due to this is more than likely to be tiny and insignificant compared to the algorithmic error made by ode45.

<RANDOM_OPINIONATED_RANT>

Fourth, WHY do you use symbolic math for this purpose?! You are doing a numerical integration, meaning, there is no analytic solution anyway. Why bother with symbolics then? Doing the integration with symbolics (through vpa even) is going to be dozens, hundreds, yes, often even thousands of times slower than keeping (or re-implementing) everything numerical (which some would argue is already slow in MATLAB compared to a bare-metal approach).

Yes, of course, for this specific, individual, isolated use case it may not matter much, but for the future I'd strongly advise you to learn to:

• use symbolics for derivations, proving theorems, simplifying expressions, ...
• use numerics to implement any algorithm or function from which actual numbers are expected.

In other words, symbolics for drafting, numerics for crunching. And exactly zero symbolics should appear in any good implementation of any algorithm.

Although it's possible to mix them to some extent, that does not mean it is a good idea to do so. In fact, that's almost never. And the few isolated cases where it is the only viable option are not a vindication of the approach. They are rare, isolated cases after all, far from the abundant norm.

For me it bears resemblance with the evil eval, with similar reasons for why it Should. Be. Avoided.

</RANDOM_OPINIONATED_RANT>

Answer 2

With the full code, it's easy to come up with a complete solution:

% Initial conditions
phi = 12.5e-3;
v0  = 300;
x0  = 0;    % (my assumption)
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 = @(t) v_T * t;
m_acc = @(t) pi * e * ro *(R_T(t)^2);
v_L = sqrt (E/ro);
R_L = @(t) v_L * t;
z   = @(t) 2 * R_L(t);
E_4 = @(t) B * ((e_r^2)* B * (0.9^(z(t)/B)-1)) /(log(0.9));

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

% Resolution of the problem
g_t = @(t)  -( phi * E * e * pi * e_r^2 * v_T / 8 + ...  % dE_1/dt
               pi * 10e9 * 2 * R_T(t) * v_T + ...        % dE_2/dt
               pi * 1e6 * e * 2 * R_T(t) * v_T );        % dE_3/dt

% The derivative of Z = [x(t); x'(t)] equals Z' = [x'(t);  x''(t)]
f = @(t,y)[y(2);
          (g_t(t) - (0.5*pi*v_T*e*ro + E_4(t)) * y(2)^2) /(y(2) * (8.33e-3 + m_acc(t)))];

% Which is readily integrated
[T,Y] = ode45(f, [0 1], [x0 v0]);

% Plot solutions
figure(1)
plot(T, Y(:,1)) 
xlabel('t [s]'), ylabel('position [m]')

figure(2)
plot(T, Y(:,2)) 
xlabel('t [s]'), ylabel('velocity [m/s]')

Results:

Note that I've not used symbolics anywhere, except to double-check my hand-derived derivatives.