Main Content

## Troubleshoot Common ODE Problems

### Error Tolerances

Question or ProblemAnswer

How do I choose the error thresholds `RelTol` and `AbsTol`?

`RelTol`, the relative accuracy tolerance, controls the number of correct digits in the computed answer. `AbsTol`, the absolute error tolerance, controls the difference between the computed answer and the true solution. At each step, the error `e` in component `i` of the solution satisfies

`|e(i)| ≤ max(RelTol*abs(y(i)),AbsTol(i))`

Roughly speaking, this means that you want `RelTol` correct digits in all solution components, excluding those smaller than the threshold `AbsTol(i)`. Even if you are not interested in a component `y(i)` when it is small, you might have to specify a value for `AbsTol(i)` that is small enough to get some correct digits in `y(i)` so that you can accurately compute more interesting components.

I want answers that are correct to the precision of the computer. Why can I not simply set `RelTol` to `eps`?

You can get close to machine precision, but not that close. The solvers do not allow `RelTol` near `eps` because they try to approximate a continuous function. At tolerances comparable to `eps`, the machine arithmetic causes all functions to look discontinuous.

How do I tell the solver that I do not care about getting an accurate answer for one of the solution components?

You can increase the absolute error tolerance `AbsTol` for this solution component. If the tolerance is bigger than the solution component, this specifies that no digits in the component need to be correct. The solver might have to get some correct digits in this component to compute other components accurately, but it generally handles this automatically.

### Problem Scale

Question or ProblemAnswer

How large a problem can I solve with the ODE suite?

The primary constraints are memory and time. At each time step, the solvers for nonstiff problems allocate vectors of length `n`, where `n` is the number of equations in the system. The solvers for stiff problems allocate vectors of length `n` but also allocate an `n`-by-`n` Jacobian matrix. For these solvers, it might be advantageous to specify the Jacobian sparsity pattern using the `JPattern` option of `odeset`.

If the problem is nonstiff, or if you are using the `JPattern` option, it might be possible to solve a problem with thousands of unknowns. In this case, however, storage of the result can be problematic. Ask the solver to evaluate the solution at specific points only, or call the solver with no output arguments and use an output function to monitor the solution.

I am solving a very large system, but only care about a few of the components of `y`. Is there any way to avoid storing all of the elements?

Yes. The `OutputFcn` option is designed specifically for this purpose. When you call the solver with no output arguments, the solver does not allocate storage to hold the entire solution history. Instead, the solver calls `OutputFcn(t,y,flag)` at each time step. To keep the history of specific elements, write an output function that stores or plots only the elements you care about.

What is the startup cost of the integration and how can I reduce it?

The biggest startup cost occurs as the solver attempts to find a step size appropriate to the scale of the problem. If you happen to know an appropriate step size, use the `InitialStep` option. For example, if you repeatedly call the integrator in an event location loop, the last step that was taken before the event is probably scaled correctly for the next integration. Type `edit ballode` to see an example.

The first step size that the integrator takes is too large, and it misses important behavior.

You can specify the first step size with the `InitialStep` option. The integrator tries this step size, then reduces it if necessary.

### Solution Components

Question or ProblemAnswer

The solution does not look like what I expected.

If your expectations are correct, then reduce the error tolerances from their default values. A smaller relative error tolerance is needed to accurately solve problems integrated over “long” intervals, as well as problems that are moderately unstable.

Check whether there are solution components that stay smaller than their absolute error tolerance for some time. If so, you are not requiring any correct digits in these components. This might be acceptable for these components, but failing to compute them accurately might degrade the accuracy of other components that depend on them.

My plots are not smooth enough.

Increase the value of `Refine` from its default of `4` in `ode45` or `1` in the other solvers. The bigger the value of `Refine`, the more output points the solver generates. Execution speed is not affected much by the value of `Refine`.

I am plotting the solution as it is computed and it looks fine, but the code gets stuck at some point.

First verify that the ODE function is smooth near the point where the code gets stuck. If it is not, then the solver must take small steps to deal with this. It might help to break the integration interval into pieces over which the ODE function is smooth.

If the function is smooth and the code is taking extremely small steps, you are probably trying to solve a stiff problem with a solver not intended for that purpose. Switch to using one of the stiff solvers `ode15s`, `ode23s`, `ode23t`, or `ode23tb`.

What if I have the final and not the initial value?

All the solvers of the ODE suite allow you to solve backward or forward in time. The syntax for the solvers is ```[t,y] = ode45(odefun,[t0 tf],y0);``` and the syntax accepts ```t0 > tf```.

My integration proceeds very slowly, using too many time steps.

First, check that `tspan` is not too long. Remember that the solver uses as many time points as necessary to produce a smooth solution. If the ODE function changes on a time scale that is very short compared to `tspan`, then the solver uses a lot of time steps. Long-time integration is a hard problem. Break `tspan` into smaller pieces.

If the ODE function does not change noticeably on the `tspan` interval, it could be that your problem is stiff. Try using one of the stiff solvers `ode15s`, `ode23s`, `ode23t`, or `ode23tb`.

Finally, make sure that the ODE function is written in an efficient way. The solvers evaluate the derivatives in the ODE function many times. The cost of numerical integration depends critically on the expense of evaluating the ODE function. Rather than recomputing complicated constant parameters at each evaluation, store them in global variables, or calculate them once and pass them to nested functions.

I know that the solution undergoes a radical change at time `t`, where `t0 <= t <= tf `, but the integrator steps past without “seeing” it.

If you know that there is a sharp change at time `t`, try breaking the `tspan` interval into two pieces, ```[t0 t]``` and `[t tf]`, and call the integrator twice or continue the integration using `odextend`.

If the differential equation has periodic coefficients or solutions, ensure the solver does not step over periods by restricting the maximum step size to the length of the period.

### Problem Type

 Can the solvers handle partial differential equations (PDEs) that have been discretized by the method of lines? Yes, because the discretization produces a system of ODEs. Depending on the discretization, you might have a form involving mass matrices, which the ODE solvers provide for. Often the system is stiff. This is to be expected if the PDE is parabolic, or when there are phenomena that happen on very different time scales such as a chemical reaction in a fluid flow. In such cases, use one of the four stiff solvers `ode15s`, `ode23s`, `ode23t`, or `ode23tb`.If there are many equations, use the `JPattern` option to specify the Jacobian sparsity pattern. This can make the difference between success and failure as it prevents the computation from being too expensive. Type ```edit burgersode``` to see an example that uses `JPattern`.If the system is not stiff, or not very stiff, then `ode23` and `ode45` are more efficient than the stiff solvers `ode15s`, `ode23s`, `ode23t`, and `ode23tb`.You can solve parabolic-elliptic partial differential equations in 1-D directly with the MATLAB® PDE solver `pdepe`. Can I integrate a set of sampled data? Not directly. Instead, represent the data as a function by interpolation or some other scheme for fitting data. The smoothness of this function is critical. A piecewise polynomial fit such as a spline can look smooth to the eye, but rough to a solver; the solver takes small steps where the derivatives of the fit have jumps. Either use a smooth function to represent the data or use one of the lower-order solvers (`ode23`, `ode23s`, `ode23t`, `ode23tb`) that is less sensitive to smoothness. See ODE with Time-Dependent Terms for an example.

Download ebook