Struggling to Solve 2nd Order ODE with Multiple Initial Values

Question

David Bloom 2023-12-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2059379-struggling-to-solve-2nd-order-ode-with-multiple-initial-values

评论： Paul 2023-12-14

在 MATLAB Online 中打开

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Sam Chak 2023-12-12

在 MATLAB Online 中打开

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 2023-12-12

编辑：Torsten 2023-12-12

@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 2023-12-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2059379-struggling-to-solve-2nd-order-ode-with-multiple-initial-values#answer_1369744

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Sam Chak 2023-12-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2059379-struggling-to-solve-2nd-order-ode-with-multiple-initial-values#answer_1370089

在 MATLAB Online 中打开

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

8 个评论
显示 6更早的评论隐藏 6更早的评论

Sam Chak 2023-12-13

在 MATLAB Online 中打开

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

Torsten 2023-12-13

编辑：Torsten 2023-12-13

在 MATLAB Online 中打开

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 2023-12-14

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.

请先登录，再进行评论。