Stiff Differential Equation solver (Euler?)

Hello, I'm looking for a fixed-step integrating function capable of solving stiff differential equations, with a really small step size (below e-14).

Is this possible? Does MATLAB have one? I couldn't make it work with the built-in ode*() functions.

Anything out there available? Maybe working on the Euler method?

NLN

[Torsten]

No. But the Euler method (although completely inadequate for stiff differential equations) is easily programmed.

[Nuno Nunes]

@Torsten

I'm not really a mastermind at either MATLAB or maths, could you help me out with it?

I understand it's a function with inputs:

• t0 (because it's in function of time for my case)
• Y0
• h
• dy

I can define h = e-16,

but what's next?

I want to work with the following obtained info:

M = [m + u_11 0 0;...
    0 m + u_22 m * x_G + u_26;...
    0 m * x_G + u_26 I_zz + u_66];
dV = M \ ...
    [X_A + X_q + X_P + m * y(2) * y(3) + m * x_G * (y(3)^2) ;...
    Y_A + Y_q - m * y(1) * y(3);...
    N_A + N_q - m * x_G * y(1) * y(3) ];
dy(1) = dV(1);                                                            % du
dy(2) = dV(2);                                                            % dv
dy(3) = dV(3);                                                            % dr

How do I write the Euler's method code for this?

[Torsten]

I read your other post in the forum.

You can easily output all variables you like from the ode function by just adding them as output arguments:

x0 = 1;
tspan = [0 1];
[T,Y] = ode45(@fun,tspan,x0);
for i = 1:numel(T)
    [~,par1(i),par2(i)] = fun(T(i),Y(i));
end
plot(par1,par2)

function [dy,par1,par2] = fun(t,y)
  par1 = t;
  par2 = y;
  dy = -y;
end

Maybe this makes your question here superfluous.

[Nuno Nunes]

I guess I must really be a noob at this. My StraightRun function is defined as:

function [dy,X,Y,N,X_AS,Y_AS,N_AS] = StraightRun(t,y) 

And I call it on my MainSimulation script as:

[T,Y] = ode45('StraightRun',t,y0);

Care to venture a guess as to why I cannot read any of the output variables besides dy?

FYI: dy is only being read on ode45 and doesnt' go through. What goes through is Y (having passed through ode45).

[Torsten]

I don't understand what you want to tell us with your comment.

Do you have problems calling the function "StraightRun" after ode45 has finished and get the additional outputs ?

In other words: Does the way I suggested above not work ?

You know that you'd better call as

[T,Y] = ode45(@StraightRun,t,y0);

than

[T,Y] = ode45('StraightRun',t,y0);

?

[Nuno Nunes]

What I meant was:

The ode45 runs ok and gives out the supposed results, but none of the other outputs on StraightRun:

function [dy,X,Y,N,X_AS,Y_AS,N_AS] = StraightRun(t,y) 

are readable (or callable) after ode45 runs for the time period t established.

(dy is a 1x7 double)

This is just bugging my mind like, how can the StraightRun function run through the solver, but leave out whatever other variables are not in 'dy' which contains the variables for the differential equation? what about everything else?

I've set them as an array 'info':

info(1) = V_A;
info(2) = beta_A;
info(3) = alpha_S;
info(4) = X;
info(5) = Y;
info(6) = N;
info(7) = X_AS;
info(8) = Y_AS;
info(9) = N_AS;

And I've even set 'info' as an output

function [dy,info] = StraightRun(t,y) 

But after the ode45 runs, when I call any of the indices in 'info', I just get that it's not defined.

Breaking it down:

I want the same thing I already am able to do to Y: plot some of its indices like:

figure % 1                                                     % Ship Trajectory
box on;
grid on;
hold on;
axis([0 2000 -700 700])
plot(Y(:,4) , Y(:,5));
xlabel('X pos [m]');
ylabel('Y pos [m]');
legend('Ship Trajectory in Straight Run');

But I want to do it to variables (in 'info') that are not transformed in ode45.

The thing is, I can't just run StraightRun (outside of ode45) for t times, because what I intend to obtain is dependent on the heading and velocity of the ship, which is something calculated in ode45.

For example

info(1) = V_A

In which V_A is the velocity of the air, which is dependent on the ship's velocity and heading and also on another set parameter on StraightRun which is the True Wind Velocity.

What I need is to read at every time step some variables working on StraightRun without having them pass through ode45 (because they are not part of a differential equation).

[Torsten]

As you can easily see from the code I included, you must call "StraightRun" for all the times returned from ODE45 with the results for t and y you got. Then your variables that you included in the output list will be recomputed for these t and y values you supplied.

[Nuno Nunes]

I see,

I've inserted this code:

            for i = 1:numel(T)
                [information] = StraightRun(T(i),Y(i));
            end

But it just returns the error

            Index exceeds the number of array elements. Index must not exceed 1.

(Renamed 'info' to 'information' to avoid problems with built-in 'info' function)

[Walter Roberson]

            for i = 1:numel(T)
                [~, information(i,:)] = StraightRun(T(i),Y(i));
            end

[Nuno Nunes]

Same thing happens @Walter Roberson

Index exceeds the number of array elements. Index must not exceed 1.

Error in MainSimulation (line 44)
        [~,information(i,:)] = StraightRun(T(i),Y(i));
 

[Torsten]

I don't know the size of the array "information" that you return from StraightRun.

Adapt it according to what you have programmed.

for i = 1:numel(T)
    [~, information] = StraightRun(T(i),Y(i,:));
end

[Nuno Nunes]

The array 'information' goes up to 9

information(1) = V_A;
information(2) = beta_A;
information(3) = alpha_S;
information(4) = X;
information(5) = Y;
information(6) = N;
information(7) = X_AS;
information(8) = Y_AS;
information(9) = N_AS;

[Torsten]

If "information" in "StraightRun" is a row vector,

for i = 1:numel(T)
    [~, information(i,:)] = StraightRun(T(i),Y(i,:));
end

Of course, your function StraightRun must have the form

function [dy,information] = StraightRun(t,y)

I hope all variables involved (V_A,beta_A,alpha_S,X,Y,N,X_AS,Y_AS,N_AS) are scalar values.

[Nuno Nunes]

I don't knwo if you're used to hearing this (pbbly are, considering that you help people), but here it goes:

@Torsten I Love you!

IT FINALLY WORKS!!

Still not sure what I'm going to do about the whole Unable to meet integration tolerances without reducing the step size below the smallest value allowed, issue, I guess I'll just tweek some of the physical parameters of the ship.

Answer 1

[John D'Errico]

You SERIOUSLY do not want to use a standard Euler's method to solve a stiff ODE. You will be wasting your time. Why do you think you want to use Euler here, when better methods are available for stiff problems?

Worse, trying to use a step size of 1e-16 is just asking for numerical problems. This will NEVER be a good idea. Period.

https://en.wikipedia.org/wiki/Euler_method

Honestly, seriously, you do NOT want to use a simple forwards Euler method here. I don't know why you think you do. But you DON'T.

Having said all of that, the backwards Euler method is an option.

https://en.wikipedia.org/wiki/Backward_Euler_method

I won't write the code for you. But the backwards (implicit) Euler method should generally be stable for stiff problems. You may still need a fine step size, but 1e-16 is just obscene.

Do some reading before you proceed, if you really think you need to write this yourself.

Having said all of that, why in the name of god and little green apples do you want to write an ODE solver code yourself? This is especially true if you don't even know the basics of these codes? Use existing code when it is available. Do you think you will write better code than that from professionals who know very well how to do the numerical analysis? Never look to write your own code, unless it is a homework assignment.

In this case, you will want to use tools like ode15s or ode23s.

help ode15s

 ODE15S Solve stiff differential equations and DAEs, variable order method.
    [TOUT,YOUT] = ODE15S(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates 
    the system of differential equations y' = f(t,y) from time T0 to TFINAL 
    with initial conditions Y0. ODEFUN is a function handle. For a scalar T 
    and a vector Y, ODEFUN(T,Y) must return a column vector corresponding 
    to f(t,y). Each row in the solution array YOUT corresponds to a time
    returned in the column vector TOUT.  To obtain solutions at specific
    times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN =
    [T0 T1 ... TFINAL].     
    
    [TOUT,YOUT] = ODE15S(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
    integration properties replaced by values in OPTIONS, an argument created
    with the ODESET function. See ODESET for details. Commonly used options
    are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector
    of absolute error tolerances 'AbsTol' (all components 1e-6 by default).  
    If certain components of the solution must be non-negative, use
    ODESET to set the 'NonNegative' property to the indices of these
    components.  The 'NonNegative' property is ignored for problems
    where there is a  mass matrix. 
    
    The Jacobian matrix df/dy is critical to reliability and efficiency. Use
    ODESET to set 'Jacobian' to a function handle FJAC if FJAC(T,Y) returns 
    the Jacobian df/dy or to the matrix df/dy if the Jacobian is constant. 
    If the 'Jacobian' option is not set (the default), df/dy is approximated 
    by finite differences. Set 'Vectorized' 'on' if the ODE function is coded 
    so that ODEFUN(T,[Y1 Y2 ...]) returns [ODEFUN(T,Y1) ODEFUN(T,Y2) ...]. 
    If df/dy is a sparse matrix, set 'JPattern' to the sparsity pattern of
    df/dy, i.e., a sparse matrix S with S(i,j) = 1 if component i of f(t,y)
    depends on component j of y, and 0 otherwise.    
 
    ODE15S can solve problems M(t,y)*y' = f(t,y) with mass matrix M(t,y). Use
    ODESET to set the 'Mass' property to a function handle MASS if MASS(T,Y) 
    returns the value of the mass matrix. If the mass matrix is constant, 
    the matrix can be used as the value of the 'Mass' option. Problems with
    state-dependent mass matrices are more difficult. If the mass matrix does
    not depend on the state variable Y and the function MASS is to be called
    with one input argument T, set 'MStateDependence' to 'none'. If the mass
    matrix depends weakly on Y, set 'MStateDependence' to 'weak' (the
    default) and otherwise, to 'strong'. In either case the function MASS is
    to be called with the two arguments (T,Y). If there are many differential
    equations, it is important to exploit sparsity: Return a sparse
    M(t,y). Either supply the sparsity pattern of df/dy using the 'JPattern'
    property or a sparse df/dy using the Jacobian property. For strongly
    state-dependent M(t,y), set 'MvPattern' to a sparse matrix S with S(i,j)
    = 1 if for any k, the (i,k) component of M(t,y) depends on component j of
    y, and 0 otherwise.    
 
    If the mass matrix is non-singular, the solution of the problem is
    straightforward. See examples FEM1ODE, FEM2ODE, BATONODE, or
    BURGERSODE. If M(t0,y0) is singular, the problem is a differential-
    algebraic equation (DAE). ODE15S solves DAEs of index 1. DAEs have
    solutions only when y0 is consistent, i.e., there is a yp0 such that
    M(t0,y0)*yp0 = f(t0,y0). Use ODESET to set 'MassSingular' to 'yes', 'no',
    or 'maybe'. The default of 'maybe' causes ODE15S to test whether M(t0,y0)
    is singular. You can provide yp0 as the value of the 'InitialSlope'
    property. The default is the zero vector. If y0 and yp0 are not
    consistent, ODE15S treats them as guesses, tries to compute consistent
    values close to the guesses, and then goes on to solve the problem. See
    examples HB1DAE or AMP1DAE.  
 
    [TOUT,YOUT,TE,YE,IE] = ODE15S(ODEFUN,TSPAN,Y0,OPTIONS) with the 'Events'
    property in OPTIONS set to a function handle EVENTS, solves as above 
    while also finding where functions of (T,Y), called event functions, 
    are zero. For each function you specify whether the integration is 
    to terminate at a zero and whether the direction of the zero crossing 
    matters. These are the three column vectors returned by EVENTS: 
    [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function: 
    VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration 
    is to terminate at a zero of this event function and 0 otherwise. 
    DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only 
    zeros where the event function is increasing, and -1 if only zeros where 
    the event function is decreasing. Output TE is a column vector of times 
    at which events occur. Rows of YE are the corresponding solutions, and 
    indices in vector IE specify which event occurred.    
    
    SOL = ODE15S(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
    used with DEVAL to evaluate the solution or its first derivative at 
    any point between T0 and TFINAL. The steps chosen by ODE15S are returned 
    in a row vector SOL.x.  For each I, the column SOL.y(:,I) contains 
    the solution at SOL.x(I). If events were detected, SOL.xe is a row vector 
    of points at which events occurred. Columns of SOL.ye are the corresponding 
    solutions, and indices in vector SOL.ie specify which event occurred. 
 
    Example
          [t,y]=ode15s(@vdp1000,[0 3000],[2 0]);   
          plot(t,y(:,1));
      solves the system y' = vdp1000(t,y), using the default relative error
      tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
      component, and plots the first component of the solution.
 
    See also ODE23S, ODE23T, ODE23TB, ODE45, ODE23, ODE113, ODE15I,
             ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
             ODEEXAMPLES, VDPODE, BRUSSODE, HB1DAE, FUNCTION_HANDLE.

    Documentation for ode15s
       doc ode15s

help ode23s

 ODE23S  Solve stiff differential equations, low order method.
    [TOUT,YOUT] = ODE23S(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates 
    the system of differential equations y' = f(t,y) from time T0 to TFINAL 
    with initial conditions Y0. ODEFUN is a function handle. For a scalar T
    and a vector Y, ODEFUN(T,Y) must return a column vector corresponding 
    to f(t,y). Each row in the solution array YOUT corresponds to a time 
    returned in the column vector TOUT.  To obtain solutions at specific 
    times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN = 
    [T0 T1 ... TFINAL]. 
    
    [TOUT,YOUT] = ODE23S(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
    integration properties replaced by values in OPTIONS, an argument created
    with the ODESET function. See ODESET for details. Commonly used options
    are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector
    of absolute error tolerances 'AbsTol' (all components 1e-6 by default).  
       
    The Jacobian matrix df/dy is critical to reliability and efficiency. Use
    ODESET to set 'Jacobian' to a function handle FJAC if FJAC(T,Y) returns 
    the Jacobian df/dy or to the matrix df/dy if the Jacobian is constant. 
    If the 'Jacobian' option is not set (the default), df/dy is approximated 
    by finite differences. Set 'Vectorized' 'on' if the ODE function is coded 
    so that ODEFUN(T,[Y1 Y2 ...]) returns [ODEFUN(T,Y1) ODEFUN(T,Y2) ...]. 
    If df/dy is a sparse matrix, set 'JPattern' to the sparsity pattern of
    df/dy, i.e., a sparse matrix S with S(i,j) = 1 if component i of f(t,y)
    depends on component j of y, and 0 otherwise.    
    
    ODE23S can solve problems M*y' = f(t,y) with a constant mass matrix M
    that is nonsingular. Use ODESET to set the 'Mass' property to the mass
    matrix. If there are many differential equations, it is important to
    exploit sparsity: Use a sparse M. Either supply the sparsity pattern of
    df/dy using the 'JPattern' property or a sparse df/dy using the Jacobian
    property. ODE15S and ODE23T can solve problems with singular mass
    matrices.     
    
    [TOUT,YOUT,TE,YE,IE] = ODE23S(ODEFUN,TSPAN,Y0,OPTIONS...) with the 'Events'
    property in OPTIONS set to a function handle EVENTS, solves as above 
    while also finding where functions of (T,Y), called event functions, 
    are zero. For each function you specify whether the integration is 
    to terminate at a zero and whether the direction of the zero crossing 
    matters. These are the three column vectors returned by EVENTS: 
    [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function: 
    VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration 
    is to terminate at a zero of this event function and 0 otherwise. 
    DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only 
    zeros where the event function is increasing, and -1 if only zeros where 
    the event function is decreasing. Output TE is a column vector of times 
    at which events occur. Rows of YE are the corresponding solutions, and 
    indices in vector IE specify which event occurred.    
 
    SOL = ODE23S(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
    used with DEVAL to evaluate the solution or its first derivative at 
    any point between T0 and TFINAL. The steps chosen by ODE23S are returned 
    in a row vector SOL.x.  For each I, the column SOL.y(:,I) contains 
    the solution at SOL.x(I). If events were detected, SOL.xe is a row vector 
    of points at which events occurred. Columns of SOL.ye are the corresponding 
    solutions, and indices in vector SOL.ie specify which event occurred. 
 
    Example
          [t,y]=ode23s(@vdp1000,[0 3000],[2 0]);   
          plot(t,y(:,1));
      solves the system y' = vdp1000(t,y), using the default relative error
      tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
      component, and plots the first component of the solution.
 
    See also ODE15S, ODE23T, ODE23TB, ODE45, ODE23, ODE113, ODE15I,
             ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
             ODEEXAMPLES, VDPODE, BRUSSODE, FUNCTION_HANDLE.

    Documentation for ode23s
       doc ode23s

[Nuno Nunes]

I get where you're coming from, and I honestly appreciate your concerns.

However, this is all inserted in a Master's dissertation where I have to simulate a ship's motion (amongst other things) under certain conditions and the ode*() functions available all give me the

Warning: Failure at t= ########. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.547474e-13) at time t.

I asked my Tutor for help, and he said I should just program my own fixed-step integrating function and that the method could even be Euler's.

So, that's that.

[Torsten]

@Nuno Nunes

If I were you, I'd check the equations you are trying to solve and the way you implemented them, but not try to write my own solver.

[John D'Errico]

(NO. The method CANNOT be a standard forward Euler's method. It will not be stable for a stiff problem. It COULD POSSIBLY be a backwards Euler though. That might depend on your skill in coding, and you knowledge of the numerical method.)

However, you are using MATLAB. It is silly to be forced to wite your own code to replace a part of MATLAB. Your advisor is insane. Tell them I said so, not that it will really matter in some cases.

Would your adisor tell you to implement addition and multiplication of numbers on your own in MATLAB, because this is your thesis? Sorry, but that is just plain dumb. Ok, don't tell your advisor that. Or do tell them. I don't care what they think of me anyway. :)

My point is still fully valid. If you are using MATLAB to solve this problem, then just use MATLAB! Using MATLAB PROPERLY to solve the problem means you need to use the correct code from MATLAB to solve it. In this case, the solution is to use the stiff ODE solvers already part of MATLAB itself. Don't write your own code! I you learn nothing more here, it is that.

Freekin' bleepin advisors.

Anyway, IF you are seriously given no choice, then this is YOUR thesis in the end. You would need to spend the mental energy to learn to write a BACKWARDS (implicit) Euler method. You can surely find pseudo-code for a backwards Euler.

[Nuno Nunes]

I gotta say... I like the way you express yourself @John D'Errico

My F'ing advisor is NOT a good person at all, he's just one of those academic guys that MUST have his name on everything. He totally hasn't helped me one bit, and I'm due to deliver this thesis at the end of the month (wish me luck).

I'm going to ignore what he's saying because it absolutely makes no sense that I should become a pro programmer, when nothing during the course has been aimed at improving my programming skills AT ALL.

Thanks for all your help and candor on this Question! Fr!

[John D'Errico]

If there never would have been a complaint had you been able to successfully use ode45, then there cannot possibly be a valid problem if you use ode15s. Both tools are essentially part of the very same suite of codes. The only important factor lies in knowing which code is the correct tool to solve the problem.

Answer 2

[Walter Roberson]

Code for fixed step solvers is at https://www.mathworks.com/matlabcentral/answers/98293-is-there-a-fixed-step-ordinary-differential-equation-ode-solver-in-matlab-8-0-r2012b#answer_107643

If this is not able to operate at a fine enough time step then you may need to alter the code to use the symbolic toolbox.

[Walter Roberson]

personally I think it likely that your equations, at least as implemented, have an unavoidable singularity.

• the equations might be wrong
• you might have coded them incorrectly
• the problem might possibly not be solvable using the techniques that you are using

I am a big fan of setting up the equations using the symbolic toolbox, which makes it much easier to follow the equations to be sure that they have been expressed correctly, especially if you use Livescript (better output format). Then use the work flow shown in the first example of odeFunction to convert the symbolic expressions for numeric solutions.