Struggling to Solve 2nd Order ODE with Multiple Initial Values

I'm currently trying to solve a 2nd order ODE with dsolve and cannot get it to properly output a solution...

Here are my ODE, initial conditions, and code...

ODE:

Initial Conditions:

: Where

My code:

clear all;

syms R(r)
eqn = diff(R,2) == R-R^3-(1/r)*diff(R);
V = odeToVectorField(eqn)

M = matlabFunction(V,'vars',{'r','Y'})

interval = [0 20];
yInit = [2.2 2.2];
ySol = ode45(M,interval,yInit);

figure(2);
tValues = linspace(0,20,100);
yValues = deval(ySol,tValues,1);
plot(tValues,yValues)

When I run this I don't get an error, but my yValues always have the first value equal to my initial condition and the rest are NaN's.

My questions:

1: How do I specify that my initial conditions are R'(0)=0 and R(large-number)=0 within the syntax of ODE45?

2: Since I assume the 1/r in my ODE is the cause of the NaN's, how do I get around this?

3: Is there a simpler way to solve this system?

Thanks in advance for your help.

-David

[Torsten]

Why worry ? The solution with the boundary values you specified is R(r) = 0 for all r - no need for computations.

[David Bloom]

How often do people solving different equations want to be told that the trivial solution is the answer? Stop being a troll.

[James Tursa]

@David Bloom

1) Torsten is not a troll, but is a very valuable poster on Answers, particularly for ODEs.

2) Your differential equations are not defined at your initial point r=0, since that 2nd term is 0/0, so what do you expect any solver to do with this? At least the R(r)=0 solution suggested by Torsten has a limit of 0 for that term as r-->0, but I'm not at all sure what that means given that the term at that point is undefined. Technically, I don't think the R(r)=0 solution qualifies as a solution because of this.

[David Bloom]

@Torsten My apologies, I'm used to interactions on Stack Exchange where the norm is to be mocked, trolled, or denigrated for asking a question.

I made a mistake in writing my question. There is no "known" value for , however it is usually found via the "shooting method" by comparing how different initial values of impact the convergence to zero of values at large values of r.

The equation in question describes the Townes beam profile which is a solution to the wave-equation in the presence of an intensity-dependant index of refraction.

I found a solution by postulating a value for , checking the convergence, and then iteratively changing the value. I also had to reduce my interval from [0 20] to [1e-33 20] to deal with the "blow-up" that occurs at r=0 due to the 1/r term.

[Sam Chak]

Hi @David Bloom

We can observe the CT scan of the phase plane diagrams for the 2nd-order ODE as r increases. These diagrams offer insights into the behavior of the system, including the equilibrium points and the solutions to the system of ODEs.

%% CT Scan of the phase plane diagrams for the 2nd-order ODE

[Torsten]

@David Bloom

The only possible boundary condition at r = 0 for your equation is either dR/dr = 0 or that the function is finite at r = 0 (if a value for R at r = 0 is not known). If you also impose the condition R(Inf) = 0, I doubt that a numerical method is able to find a solution apart from R = 0. Try with "bvp4c" which is the boundary value solver in the MATLAB software.

Answer 1

[Walter Roberson]

syms R(r)
eqn = diff(R,2) == R-R^3-(1/r)*diff(R);
[eqs,vars] = reduceDifferentialOrder(eqn,R(r))

eqs = 

vars = 

[M,F] = massMatrixForm(eqs,vars)

M = 

F = 

f = M\F

f = 

odefun = odeFunction(f,vars)

odefun = function_handle with value:
    @(r,in2)[in2(2,:);-(in2(2,:)-in2(1,:).*r+in2(1,:).^3.*r)./r]

interval = [0 20];
yInit = [2.2 2.2];
[tValues, yValues] = ode45(odefun,interval,yInit);
[min(tValues), max(tValues)]

ans = 1×2

0 20

[min(yValues(:,1)), max(yValues(:,1))]

ans = 1×2

2.2000 2.2000

nnz(isnan(yValues(:,1)))

ans = 1888

plot(tValues,yValues(:,1))

limit(lhs(eqn), r, 0, 'left') == limit(rhs(eqn), r, 0, 'left')

ans = 

limit(lhs(eqn), r, 0, 'right') == limit(rhs(eqn), r, 0, 'right')

ans = 

No output because the equation involves division by r and r starts out as 0.

Numeric processing of 0/0 is not at all forgiving of the fact that there might potentially be defined value there due to limit reasons. As you can see, MATLAB cannot inherently reason about the limits for those functions anyhow -- more complex analysis would have to be done to determine whether anything can be figured out about the limit of the derivative at 0

1: How do I specify that my initial conditions are R'(0)=0 and R(large-number)=0 within the syntax of ODE45?

There is no explicit syntax for that for ode45.

In terms of the processing above you would have to analyze the vars returned by reduceDifferentialOrder in order to determine which index corresponded to which initial condition.

You might also be interested in the new (R2023b) facility ode which permits a different way of constructing ODEs.

One of the features of the new facility is that it permits initial conditions to be specified and the associated times, but to solve the ode over a different time range. So for example you might specify r(0) value but ask to run the ode over time -10 to +10 (example), and it (somehow) knows how to do that without you having to call the ode solvers twice, once to work backwards from 0 to -10 and the other to work forwards from 0 to 10.

Answer 2

[Sam Chak]

Hi @David Bloom

I'm not an expert in self-trapped optical spatial solitons, but I'm attempting to reproduce the graph shown in Figure 1 in Townes et al.'s original paper: https://www.researchgate.net/publication/220030032_Self-Trapping_of_Optical_Beams.

I have also provided a set of initial values for that satisfy the boundary condition at .

%% Case 1: Reproduce the graph in Townes's paper
rspan   = linspace(1e-6, 5, 50001);
R0      = [2.2075; 0];
[r, R]  = ode45(@odefcn, rspan, R0);
figure(1)
plot(r, R(:,1)), grid on
xlabel('r'), ylabel('R(r)')

%% Case 2: Initial values of R that satisfy R(15) ≈ 0
Rinit   = [5.4637, 6.4430, 7.6445, 8.9731, 10.4083, 11.9213];
figure(2)
for j = 1:numel(Rinit)
    rspan   = linspace(1e-6, 15, 150001);
    R0      = [Rinit(j); 0];
    [r, R]  = ode45(@odefcn, rspan, R0);
    plot(r, R(:,1)), hold on
end
hold off, grid on
xlabel('r'), ylabel('R(r)')

%% Dynamics of Spatial Soliton
function dRdr = odefcn(r, R)
    dRdr    = zeros(2, 1);
    dRdr(1) = R(2);
    dRdr(2) = R(1) - R(1)^3 - (1/r)*R(2);
end

[Sam Chak]

Hi @David Bloom

It should be noted that is not a stable equilibrium point. If you simulate the dynamics for a sufficient duration, eventually, four trajectories will circle around one of the two stable equilibrium points, . Another two will circle around . The results are expected and consistent with the phase plane diagram above as r increases.

%% Case 2: Initial values of R that satisfy R(15) ≈ 0
Rinit   = [5.4637, 6.4430, 7.6445, 8.9731, 10.4083, 11.9213];
figure(2)
for j = 1:numel(Rinit)
    rspan   = linspace(1e-6, 150, 150001);
    R0      = [Rinit(j); 0];
    [r, R]  = ode45(@odefcn, rspan, R0);
    plot(R(:,1), R(:,2)), hold on
end
hold off, grid on
xlabel('R'), ylabel('dR/dr')
axis([-pi/2 pi/2 -pi/2 pi/2])

%% Dynamics of Spatial Soliton
function dRdr = odefcn(r, R)
    dRdr    = zeros(2, 1);
    dRdr(1) = R(2);
    dRdr(2) = R(1) - R(1)^3 - (1/r)*R(2);
end

[Torsten]

@Sam Chak

What is your conclusion ? Is there another condition on the solution that comes from the physical background necessary to sort out the "correct" one ? Does perhaps dR/dr < 0 or R > 0 for all r suffice ?

[Sam Chak]

Hi @Torsten

My conclusion is that the system cannot converge to as r approaches ∞. Regardless of the OP's desire, mathematical analysis from the phase plane diagram indicates that is not a steady-state. I'm not an expert, but perhaps if the OP injects energy into the optical beam system to force the system to converge to R = 0, then it is 100% achievable.

The OP probably found the ODE in this paper. However, as I'm not an optical physicist, I don't fully comprehend how the author evaluates the limit. There might be hidden information that the OP hasn't fully informed us about.

Self-Similar Optical Wave Collapse: Observation of the Townes Profile

[Torsten]

So you also think that R = 0 for all r is the only possible solution ?

[Sam Chak]

Hi @Torsten

Yes, I 100% agree that the trivial solution is the only possible solution. The OP can test this and observe the results. Personally, I prefer to conduct analyses, which is why I initially posted the phase plane diagrams (as r increases) for the OP to visualize the possible trajectory solutions.

rspan   = linspace(1e-6, 100, 100001);
R0      = [0; 0];
[r, R]  = ode45(@odefcn, rspan, R0);

figure(1)
subplot(211)
plot(r, R(:,1)), grid on
xlabel('r'), ylabel('R(r)')
subplot(212)
plot(R(:,1), R(:,2), '.', 'markersize', 12), grid on
xlabel('R'), ylabel('dR/dr')

However, I also believe that the OP is attempting to reproduce the results as guided in the optical beam paper. At one point, I wondered if there was a typo in the sign of the differential equation. I genuinely wanted to help, but this is as far as I can go.

Phase plane diagrams for

[x, y] = meshgrid(-pi/2:0.05:pi/2);
r      = [0.1, 1, 10];
for i = 1:numel(r)
    u     =   y;
    v     = - x - x.^3 - (1/r(i))*y;
    figure(i+1)
    l     = streamslice(x, y, u, v);
    xlabel('R'), ylabel('R^{\prime}')
    title(['r = ', num2str(r(i))])
    axis tight
end

%% Dynamics of Spatial Soliton
function dRdr = odefcn(r, R)
    dRdr    = zeros(2, 1);
    dRdr(1) = R(2);
    dRdr(2) = R(1) - R(1)^3 - (1/r)*R(2);
end

[Paul]

Can someone clarify how that boundary condition at r = 0 works?

Suppose I had a function R(r) that I claimed is a solution to the differential equation. The next step would be to show that R(r) satisfies the differential equation and satisfies the boundary conditions. In order to show the latter, is it true that not only

(1) dR(r)/dr | r=0 = 0

but also that

(2) dR(r)/dr has to go to zero at least as fast as r?

That is, the boundary condition at r = 0 has to be viewed as a limit?

Is (2) actually an implicit boundary condition on the solution or am I thinking about the problem incorrectly?

If (2) is a condition, does the form of the diferential equation automatically ensure that (2) will be satisfied for any solution that actually pushes the boundary condition on dR(r)/dt to r = 0+ (like @Sam Chak did by pushing it to r = 1e-6)?

[Torsten]

I'm also not 100% sure.

I think the condition du/dr = 0 | r=0 for equations with a differential term of the form

1/r^m * d/dr (r^m * du/dr) (m=1 for a cylindrical coordinate system, m=2 for a radial coordinate system)

is kind of an artificial condition to sort out solutions that are unbounded at r = 0.

Consider e.g.

syms r u(r)
eqn = 1/r * diff(r*diff(u,r)) == 0;

with solution

sol = dsolve(eqn)

sol = 

dsol = diff(sol)

dsol = 

A bounded solution must have C_2 = 0. To ensure this, we demand dsol = 0 at r= 0 which forces C_2 to be 0.

The condition is also physically reasonable since du/dr = 0 at r = 0 means: Nothing of the quantity u should be able to vanish over the axis r = 0. And this makes sense since the axis - seen as a cylindrical area - has area 0.

[Paul]

This example disproves my hypothesized condition (2) and is making me inclined to think that the boundary condition should be interpreted as a condition at r = 0+. But I'm still not sure.