Example 6.1.2

The equation $u'=\sin[(u+t)^2]$ also has a solution that can be found numerically with ease, even though no formula exists for its solution.

The solve_ivp function from scipy.integrate package offers solvers for a variety of initial-value problems. Note that even though this is a problem for a scalar function $u(t)$, we had to set the initial condition as a "one-dimensional vector."

You can see above that the solution was not computed at enough points to make a smooth graph. There is a way to request output at times of your choosing.

Example 6.1.3

We return to the equation $u'=\sin[(u+t)^2]$ to inspect the output a bit more closely. With a different calling sequence, we get a different form of the solution.

The t and y fields are as before. But there is another field that can be called like any other function of $t$.

As you might guess, the solution object performs some interpolation on some discrete solution values This chapter is about how the discrete $t$ and $u$ values are computed.

Example 6.1.5

The equation $u'=(u+t)^2$ gives us some trouble.

We requested the solution for $0\le t \le 1$, but the integration did not complete that job.

A plot shows that the solution appears to have blown up along a vertical asymptote. This can happen in nonlinear ODEs.

Example 6.1.6

Consider the ODEs $u'=u$ and $u'=-u$. In each case we compute $\partial f/\partial u = \pm 1$, so the condition number bound is $e^{(b-a)}$ in both problems. However, they behave quite differently. In the case of exponential growth, $u'=u$, the bound is the actual condition number.

But with $u'=-u$, solutions actually get closer together with time.

In this case the actual condition number is one, due to the difference of solutions at the initial time.

Example 6.2.1

We consider the IVP $u'=\sin[(u+t)^2]$ over $0\le t \le 4$, with $u(0)=-1$.

We could define a different interpolant to get a smoother picture above, but the derivation assumed the piecewise linear interpolant, so it is the most meaningful one. We can instead request more steps to make the interpolant look smoother.

Increasing $n$ changed the solution noticeably. Since we know that interpolants and finite differences become more accurate as $h\to 0$, we should expect that from Euler's method too.

We don't have an exact solution to compare to, so we will use solve_ivp to construct an accurate solution.

Now we can perform a convergence study.

The error is almost perfectly halved at each step, so we expect that a log-log plot will reveal first-order convergence.

Example 6.3.2

Example 6.3.3

We encode the predator–prey equations via a function.

Note that the function accepts four inputs. In order to use it, we also have to assign values to alpha and beta.

To solve the IVP we must also provide the initial condition, which is a 2-vector here, and the interval for the independent variable.

The plot above lacks smoothness because it uses just the discrete values, without interpolating between them.

When there are just two components, it's common to plot the solution in the phase plane, i.e., with $u_1$ and $u_2$ along the axes and time as a parameterization of the curve.

From this plot we can deduce that the solution approaches a periodic one, which in the phase plane is reprepresented by a closed loop.

Example 6.4.1

We consider the IVP $u'=\sin[(u+t)^2]$ over $0\le t \le 4$, with $u(0)=-1$.

We use solve_ivp to construct an accurate approximation to the exact solution.

Now we perform a convergence study of our two Runge--Kutta implementations.

The amount of computational work at each time step is assumed to be proportional to the number of stages. Let's compare on an apples-to-apples basis by using the number of $f$-evaluations on the horizontal axis.

The fourth-order variant is more efficient in this problem over a wide range of accuracy.

Example 6.5.1

Let's run adaptive RK on $u'=e^{t-u\sin u}$.

The solution makes a very abrupt change near $t=2.4$. The resulting time steps vary over three orders of magnitude.

If we had to run with a uniform step size to get this accuracy, it would be

On the other hand, the average step size that was actually taken was

We took fewer steps by a factor of almost 1000! Even accounting for the extra stage per step and the occasional rejected step, the savings are clear.

Example 6.5.2

Example 6.7.1

We consider the IVP $u'=\sin[(u+t)^2]$ over $0\le t \le 4$, with $u(0)=-1$.

We use a solve_ivp to construct an accurate approximation to the exact solution.

Now we perform a convergence study of the AB4 code.

The method should converge as $O(h^4)$, so a log-log scale is appropriate for the errors.

Example 6.7.2

The following simple ODE uncovers a surprise.

We will solve the problem first with the implicit AM2 method using $n=200$ steps.

Now we repeat the process using the explicit AB4 method.

Once the solution starts to take off, the AB4 result goes catastrophically wrong.

We hope that AB4 will converge in the limit $h\to 0$, so let's try using more steps.

So AB4, which is supposed to be more accurate than AM2, actually needs something like 8 times as many steps to get a reasonable-looking answer!

Example 6.8.1

Consider the ridiculously simple IVP $u'=u$, $u(0)=1$, whose solution is $e^t$.

Let's apply the LIAF method to this problem for varying fixed step sizes. We'll measure the error at the time $t=1$.

The error starts out promisingly, but things explode from there. A graph of the last numerical attempt yields a clue.

It's clear that the solution is growing exponentially in time.