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

13 次查看（过去 30 天）
Chris Hooper
Chris Hooper 2025-4-4
编辑： James Tursa 2025-4-5
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')
  5 个评论
显示 3更早的评论
Torsten
Torsten 2025-4-5
编辑：Torsten 2025-4-5
@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
James Tursa 2025-4-5
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2025-4-4
编辑：Torsten 2025-4-4
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.
  4 个评论
显示 2更早的评论
Torsten
Torsten 2025-4-5
编辑：Torsten 2025-4-5
@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 ).

请先登录，再进行评论。

更多回答（1 个）

Sam Chak
Sam Chak 2025-4-5
Ran in:
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)
  4 个评论
显示 2更早的评论
Sam Chak
Sam Chak 2025-4-5
编辑：Sam Chak 2025-4-5
Ran in:
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')

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by