ode45 not able to solve one particular IVP

Hello there,

the ode (attached) can be solved at three initial conditions (i=3 in the code). However when I change the i to i=4 the code fails. Is there anyway ode45 can be augmented with options to get the solution for the fourth initial condition. Other numerical programs like ode23, etc do not work either. I tried to use dsolve to get the symbolic solution, but output function involves bessel equation so that the plots look very weired. The problem is not a misprint because I am able to use RungeKutta4 to get the fourth solution on positive x-axis but the solution on the negative x-axis is still missing (trying to fix it).

please suggest how can the ode be resolved using ode45/23 etc.!

Thanks.

-Saurav Parmar

% © S. Parmar
% last modified 12-12-2023
% Zill Problem 2 page 38
clear
clc


y = -3:0.4:3;                    % pick values of y
x = -3:0.4:3;                    % pick values of x
[xg,yg] = meshgrid(x,y);         % create a 2-d array of samples of x and y
dy = xg.^2 - yg.^2;              % create dy/dx
dx = ones(size(dy));             % create corresponding time coordinates

figure
ax = gca;
quiver(xg,yg,dx,dy,1.0)                       % quiver displays the directional arrows
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
xlabel('$$\rm x $$','Interpreter','latex')
ylabel('$$ \rm y $$','Interpreter','latex')
axis tight
ax.FontSize = 16.5;
set(gca,'FontName','Times')
title('$$ \rm dy/dt=x^2 - y^2 $$','Interpreter','latex')

hold on

clear x y

syms y(x) C

% part (b)

xf = [-3.1 3.1];
xi = [-2 3 0 0];
yi = [1 0 0 2];



for i = 1:4
  for j = 1:2
      [xs,ys] = ode23(@ODEexpn,[xi(i) xf(j)],yi(i));       % get the symbolic solution
      disp(xs)
      disp(ys)
      plot(xs,ys,'color','k',LineWidth=2.5)                % plot the solution
      % hold on
  end
end



fig = gcf;
fig.Units = 'centimeters';
fig.Position = [1 2 25 15.4];
set(gcf,'PaperUnits','centimeters','PaperPosition',[1 2 25 15.4])


function dydx = ODEexpn(x,y)
dydx = x^2 - y^2;
end

[Steven Lord]

What does "fails" mean in this context?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.

[Saurav]

I get something which is not expected. I get three solutions as in ode45.png with ODE45 (and that is fine) and the 4th solution comes up with RK4 as shown in ode45+rk4.png. however when I use ODE45 for the 4th solution, I get ode45_error.png. clearly, there something going wrong in ODE45.

image(imread('ode45.png'))

Error using imread>get_full_filename (line 633)
Unable to find file "ode45.png".

Error in imread (line 395)
fullname = get_full_filename(filename);

[Torsten]

@Rena Berman

It's not possible to reference attached files in Answer (see above). Are you still working on the problem ?

[Saurav]

yes the problem is still there; it is not fixed yet. if I run the original code with i = 1:3, I get three solution as shown in ode45.png. but if i run the code with i = 1:4 I get the result as shown in ode45_errror.png. The png files are not being read in the original code itself.

[Torsten]

Could you include the png-files as images using the "image" button in the menu bar and not as files ?

[Paul]

The .png files are viewable if you click on them as a link, if that helps. At least they are for me.

[Torsten]

@Paul

Yes, that works. Thank you.

@Saurav

Seems your differential equation cannot be solved in [0 -3.1] if you impose y(0)=2 as initial condition.

syms x y(x)
eqn = diff(y,x) == x^2-y^2;
cond = y(0)==2;
sol = dsolve(eqn,cond);
fplot([real(sol),imag(sol)],[-3.1 3.1])

[Saurav]

yep I understand. But why am I not able to get the solution on [0 3.1] using ode45? Did you try the dsolve on first three conditions?..I am sure you will not be able to solve two of them using dsolve (you may try and see for yourself). As I understand this problem, we cannot use dsolve to solve this DE; ultimately numerical package has to be used.

[Torsten]

But why am I not able to get the solution on [0 3.1] using ode45?

But it works on [0 3.1]:

ode =@(x,y) x^2-y^2;
tspan = [0 3.1];
y0 = 2;
[T,Y] = ode45(ode,tspan,y0);
plot(T,Y)

The only constellation that it doesn't work is on [0 -3.1].

I am sure you will not be able to solve two of them using dsolve (you may try and see for yourself). As I understand this problem, we cannot use dsolve to solve this DE; ultimately numerical package has to be used.

"dsolve" gives two solutions - maybe that's why the difficulties arise:

syms x y(x)
eqn = diff(y,x) == x^2-y^2;
sol = dsolve(eqn);
char(sol)

ans = '[(besselj(1/4, (x^2*1i)/2)/(2*x^(1/2)) + x^(1/2)*(x*besselj(-3/4, (x^2*1i)/2)*1i - besselj(1/4, (x^2*1i)/2)/(2*x)) + (C1*bessely(1/4, (x^2*1i)/2))/(2*x^(1/2)) + C1*x^(1/2)*(x*bessely(-3/4, (x^2*1i)/2)*1i - bessely(1/4, (x^2*1i)/2)/(2*x)))/(x^(1/2)*besselj(1/4, (x^2*1i)/2) + C1*x^(1/2)*bessely(1/4, (x^2*1i)/2)); (bessely(1/4, (x^2*1i)/2)/(2*x^(1/2)) + x^(1/2)*(x*bessely(-3/4, (x^2*1i)/2)*1i - bessely(1/4, (x^2*1i)/2)/(2*x)))/(x^(1/2)*bessely(1/4, (x^2*1i)/2))]'

[Saurav]

yep...the bessel equation complicates the matter and possibly the nature of solution is different on [-3.1 0].

[Rena Berman]

(Answers dev) @Torsten, I have reported the problem internally and we are still working on it.

[Tushal Desai]

(Answers dev) The uploaded files should now be accessible when running your code. The issue has been fixed on our end. Thanks for reporting.

[Torsten]

Thank you for your efforts.

Answer 1

[Alan Stevens]

How about:

y = -3:0.4:3;                    % pick values of y
x = -3:0.4:3;                    % pick values of x
[xg,yg] = meshgrid(x,y);         % create a 2-d array of samples of x and y
dy = xg.^2 - yg.^2;              % create dy/dx
dx = ones(size(dy));             % create corresponding time coordinates

figure
ax = gca;
quiver(xg,yg,dx,dy,1.0)                       % quiver displays the directional arrows
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
xlabel('$$\rm x $$','Interpreter','latex')
ylabel('$$ \rm y $$','Interpreter','latex')
axis tight
ax.FontSize = 16.5;
set(gca,'FontName','Times')
title('$$ \rm dy/dt=x^2 - y^2 $$','Interpreter','latex')

hold on

% part (b)

xf = [3.1 -3.1];
xi = [-2 3 0 0];
yi = [1 0 0 2];


for i = 1:4
    y0 = yi(i);
  for j = 1:2
      if i==4 && j==2, y0 = -2; end
      [xs,ys] = ode45(@ODEexpn,[xi(i) xf(j)],y0);       
      plot(xs,ys,'color','k','LineWidth',2.5)        % plot the solution
  
  end
end


fig = gcf;
fig.Units = 'centimeters';
fig.Position = [1 2 25 15.4];
set(gcf,'PaperUnits','centimeters','PaperPosition',[1 2 25 15.4])
axis([-3.1 3.1 -3 3])

function dydx = ODEexpn(x,y)
dydx = x^2 - y^2;
end

[Saurav]

This is very insightful. so it looks like ode45 is working. but the nature of solution is different from what i expected. but thanks a lot. it looks like a starting point for the 4th initial condition.