how to deal with discontinuities in solutions to differential equations using ode45

Motivation: I am trying to learn how to use ode45 for fun after studying some math (because it seems engineering eventually makes you forget math and because I didn't have access to matlab in college many years ago).

Question: Given a differential equation with a discontinuity such as y'=y^2 at initial condition y(0)=3, how can I use ODE45 to approximate the solution before and after the discontinuity? I am able to get the right answer by restarting the integration at the discontunity with a new initial condition, but this requires me to know what the discontinuity looks (is it +inf or -inf taking the right hand limit)? Is there a way to use ODE45 that does not require this knowledge of the exact solution?

%The numerical solution:
f=@(t,y) y.^2;
solution=ode45(f,[0,1/3],3);
plot(solution.x,solution.y); hold on
solution2=ode45(f,[1/3,10],-1000000); %choose a number near -inf because we know the exact solution
plot(solution2.x,solution2.y); hold on

% plot(t,y_ns_r);
ylabel('y(t)')
xlabel('t')
ylim([-10 10])
%The exact solution:
t=[t_l;transpose(t_l(end)+eps:0.1:10)];
y_es=1./(-t+1/3);

plot(t,y_es,'o')

[Walter Roberson]

syms y(t)
dy = diff(y,t);
eqn = dy == y^2;
ic = y(0) == 3;
sol = dsolve(eqn, ic)

sol = 

%cross check
subs(eqn, y, sol)

ans(t) = 

There does not appear to be a discontinuity for this equation and this initial condition.

It also appear fine for other initial conditions

syms y0
ic2 = y(0) == y0;
sol2 = dsolve(eqn, ic2)

sol2 = 

limit(sol2, y0, 0)

ans = 

0

limit(sol2, y0, -inf)

ans = 

limit(sol2, y0, inf)

ans = 

[Chris Hooper]

Hi Walter, thanks for the response. I must be misunderstanding somehow. This is the graph and error message that I get:

Is the solution not discontinous at t=1/3?

[James Tursa]

@Chris Hooper Yes, it is discontinuous at t = 1/3.

[Torsten]

@Chris Hooper Strictly speaking, your solution has a pole (singularity) at t = 1/3 (and thus is not defined at this point). When talking about "continuous" or "discontinuous" at a point, the function must be defined there and have a value. This is not the case here.

[James Tursa]

See my comments below about what a "solution" is, but for clarity, by "discontinuous" in my comment above I was only referring to the function itself aside from the "solution" issue.

Answer 1

[Torsten]

Usually, the differential equations to be solved are much more complicated than yours which makes it impossible to foresee possible "discontinuities" of the solution. So the only way to proceed is to integrate and analyze the results up to the value of the independent variable where the solver gets stuck. And it's usually undecidable whether it's a real problem inherent to the equations or a weakness of the numerical method in use.

[Chris Hooper]

Thanks for the answer! Is it also undecidable whether to restart the integration from positive or negative infinity (a large number in stand in for ininity) after a suspected discontinuity?

[James Tursa]

@Chris Hooper Technically, there is no single solution that crosses the t = 1/3 boundary. You would be working with two separate solutions to the differential equation with two different initial (or limiting) conditions. If that makes sense for your application, so be it. But the solution that starts with y(0) = 3 never crosses t = 1/3, period. That solution is only defined on the interval t ∈ [0,1/3)

[Torsten]

@Chris Hooper You must know a new initial value to the right of the singularity in order to restart the integration. The two branches are not connected - so it would be like a miracle if you could catch the "correct" branch to continue the solution left of the singularity ( e.g. by choosing (2/3,-3) as new initial starting point in your case ).

[Chris Hooper]

Thanks for the answers, I think I've got it.

Answer 2

[Sam Chak]

Hi @Chris Hooper

To answer your question posed in the title, when simulating systems using ode45() or any other numerical solvers, it is generally not possible to visually determine whether the solutions have discontinuities based solely on the graphs. Here are two examples in which both solutions appear to stop at second, despite we tell the solver to integrate from 0 to 1 second.

However, in the first case, the solution actually exhibits no discontinuity, whereas in the second case (which is identical to yours), it does. If @Walter Roberson had not provided the analytical solution for the second case using dsolve(), we would not have been aware that the solution has a discontinuity at second.

Fortunately, the ode45 solver includes a warning system to alert the user if it is unable to continue with the numerical integration without reducing the step size below the minimum allowable value.

Case #1: Unstable 1st-order system,

%% Case #1: Unstable 1st-order system
[t, y]  = ode45(@(t, y) 2097*y, [0, 1], 3);    % sim from 0 to 1 second
idx     = numel(y(~isnan(y)));
tFinal  = t(idx)

tFinal = 0.3333

figure
plot(t, y), grid on, xlabel('Time / s')
xline(1/3, '--r');
title({'Solution of $\frac{dy}{dt} = 2097 y,\; y(0) = 3$'}, 'interpreter', 'latex', 'fontsize', 16)

●

Case #2: Semi-stable 1st-order system,

%% Case #2: Semi-stable 1st-order system
[t, y]  = ode45(@(t, y) y^2, [0, 1], 3);    % sim from 0 to 1 second

Warning: Failure at t=3.333305e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.881784e-16) at time t.

figure
plot(t, y), grid on, xlabel('Time / s')
xline(1/3, '--r');
title({'Solution of $\frac{dy}{dt} = y^{2},\; y(0) = 3$'}, 'interpreter', 'latex', 'fontsize', 16)

[Sam Chak]

Hi @Chris Hooper

From the figure below, it is evident that there are infinitely many solutions for the 1st-order differential equation, . However, there are three patterns of solutions, depending on the initial value, .

• If , the solution will grow unbounded. (red curve)
• If , the solution will remain at 0 indefinitely. (amber line)
• If , the solution will converge asymptotically to the asymptote , but will never reach exactly 0 in finite time. (green curve)

[Chris Hooper]

@Sam Chak Great examples thank you! Now I am excited about building plots with families of solutions. How were you able to create that one?

[Sam Chak]

Hi @Chris Hooper

Please take another look at the updated family of solution curves. The red curve and the green curve are clearly not connected. However, they coincidentally share the same analytical expression: . Thus, each positive branch has its negative counterpart, or vice versa. In fact, they all belong to the same family of solutions, expressed as, , where is the initial value at time, .

%% Green curve
y0  = 3;                % initial value at time t = 0
t   = 0.4;
y   = y0/(1 - y0*t)     % initial value to restart after discontinuity t = 1/y0 at t = 0.4

y = -15.0000

Using for-loop approach to generate the solution branches of

t0  = linspace(0, 0.7, 8);      % can use random numbers or integers
y0  = linspace(-15, 15, 31);    % can use random numbers or integers

hold on
for i = 1:numel(y0)
    for j = 1:numel(t0)
        [t, y]  = ode45(@(t, y) y^2, [t0(j), 1], y0(i));
        warning('off');
        plot(t, y, 'color', "#0072BD");
    end
end
hold off
axis tight
axis([0.0, 0.8, -15, 15])
xlabel('t'), ylabel('y(t)')
title('Family of solution curves')

[Chris Hooper]

Ah thank you!