Stiff Differential Equation solver (Euler?)

28 views (last 30 days)
Nuno Nunes
Nuno Nunes on 6 Oct 2022
Commented: Nuno Nunes on 8 Oct 2022
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
  14 Comments
Nuno Nunes
Nuno Nunes on 8 Oct 2022
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.

Sign in to comment.

Answers (2)

John D'Errico
John D'Errico on 7 Oct 2022
Edited: John D'Errico on 7 Oct 2022
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.
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.
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
  5 Comments
John D'Errico
John D'Errico on 7 Oct 2022
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.

Sign in to comment.


Walter Roberson
Walter Roberson on 8 Oct 2022
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.
  1 Comment
Walter Roberson
Walter Roberson on 8 Oct 2022
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!