How to deal with zero in NDSolve in mathematica?

I would like to solve the following ODEs $$\begin{cases} x'(t)&=y\\ y'(t)&=-y(t)/t-e^{x(t)},\\ x(0)&=1,\\y(0)&=0, \end{cases}$$

(EDIT :

The second equation used to be $y'(t) = -x(t)/t - e^{x(t)}$).

with the following code:

NDSolve[{
{x'[t] == y[t], y'[t] == -y[t]/t - Exp[x[t]]},
{x == 1, y == 0}}
,
{x, y}, {t, 0, 1}]

I get the warning:

Power::infy: "Infinite expression 1/0. encountered."

How do I fix it?

Note that, in some cases, a "small" change of the initial values will make a "large" change in in the solution, so substituting $0$ with a very small number, say $10^{-100}$ is not "good" enough.

Are there other ways to fix this?

• Maybe you could consider a variable substitution? – J. M. will be back soon May 14 '13 at 16:05
• are you sure this is a sound set of initial conditions? – Spawn1701D May 14 '13 at 16:19
• @whuber, maybe $\lim\limits_{t\to 0} t y^\prime(t)$ is an indeterminate form? – J. M. will be back soon May 14 '13 at 16:43
• Probably the initial condition is a limiting initial condition, this is common when the general solution is expressed by transcendental and special functions. Perhaps you should perturb t by $\epsilon$ and then take $\epsilon\rightarrow 0$. – Spawn1701D May 14 '13 at 17:21
• That completely changes the question! – whuber May 15 '13 at 12:40

I think your initial condition is singular. In order to solve the ODEs, consider $(x(t), y(t))$ as a planar curve, we may try changing the parameter $t$ to the arc length parameter $s$:

$$\left\{\begin{split} \frac{\mathrm{d}x}{\mathrm{d}s}=\frac{x'(t)}{\sqrt{x'(t)^2+y'(t)^2}}\\ \frac{\mathrm{d}y}{\mathrm{d}s}=\frac{y'(t)}{\sqrt{x'(t)^2+y'(t)^2}}\\ \frac{\mathrm{d}t}{\mathrm{d}s}=\frac{1}{\sqrt{x'(t)^2+y'(t)^2}} \end{split}\right.$$

Without loss of generality, we can set $t(s=0)=0$. And assume $t(s)=0$ only for countable $s$, so instead of the original ODEs, we now have:

$$\left\{\begin{split} x'(s)&=\frac{y(s)t(s)}{\sqrt{\left[\mathrm{e}^{x(s)}t(s)+x(s)\right]^2+t(s)^2 y(s)^2}}\\ y'(s)&=\frac{\mathrm{e}^{x(s)}t(s)+x(s)}{\sqrt{\left[\mathrm{e}^{x(s)}t(s)+x(s)\right]^2+t(s)^2 y(s)^2}}\\ t'(s)&=\frac{t(s)}{\sqrt{\left[\mathrm{e}^{x(s)}t(s)+x(s)\right]^2+t(s)^2 y(s)^2}}\\ x(0)&=1\\ y(0)&=0\\ t(0)&=0\\ \end{split}\right.$$

Solve it in Mathematica:

sol = NDSolve[{
Derivative[x][s] == (y[s] t[s])/Sqrt[(E^x[s] t[s] + x[s])^2 + t[s]^2 y[s]^2],
Derivative[y][s] == -((E^x[s] t[s] + x[s])/Sqrt[(E^x[s] t[s] + x[s])^2 + t[s]^2 y[s]^2]),
Derivative[t][s] == t[s]/Sqrt[(E^x[s] t[s] + x[s])^2 + t[s]^2 y[s]^2],
x == 1, y == 0, t == 0},
{x, y, t}, {s, 0, 1}] Plot the sol:

ParametricPlot[Evaluate[{{s, x[s]}, {s, y[s]}, {s, t[s]}} /. sol[]],
{s, 0, 1},
PlotRange -> All, AspectRatio -> 1/GoldenRatio,
PlotStyle -> {Blue, Darker[Green], Directive[Red, Dashed]},
Frame -> True, Axes -> False] We can see the singularity is shown clearly. And $t$ is constantly equals to $0$ on whole $s\in[0,1]$ interval, which is inconsistent with our previous assumption. So the original system might be ill-defined.

• The reparametrize the solution curve in arc length is a good idea! – van abel May 15 '13 at 3:31
• @vanabel As long as your curve is regular, it can always be parametrized by the arc length. However your original system seems not lead to a regular solution curve, so I guess your system is ill-defined. – Silvia May 15 '13 at 3:41

- update -

@whuber gives insightful comment and I agree I should mention here the following.

Below we consider a perturbed form (solvable exactly in Bessel functions) of original equation up to linear term. This allows us to understand behavior of the system around t ~ 0.

- original -

The following shows that under some general assumptions there is no solution to the system. I am not sure though if these assumptions are those that user imposes.

When t ~ 0 you claim we should have x ~ 1. If this is true we can expand Exp[x[t]] around x~1 or s~0 if x = 1+s. The resulting equation should be close to your solution (if it exists) at t ~ 0 (note I excluded y completely):

sol[t_] = First[s[t] /. DSolve[{s''[t] == -(1 + s[t])/t, s == 0}, s, t]]
-1 - Pi Sqrt[t] BesselY[1, 2 Sqrt[t]] + Sqrt[t] BesselJ[1, 2 Sqrt[t]] C

Now you have to apply also s' == 0 which is easy to prove cannot be held:

Limit[D[sol[t], t], t -> 0]
Infinity

You could also see it in a simpler way as Mathematica reterns empty solution set:

DSolve[{s''[t] == -(1 + s[t])/t, s == 0, s' == 0}, s, t]
{}
• +1 This argument can be made rigorous by considering the approximate differential equation as a perturbation of the original. Applying DSolve with initial conditions is not very persuasive. It's more insightful to omit initial conditions: the solution will include a Bessel function of the second kind. This captures the singular behavior at the origin and approximates the singular behavior of the original equation there, too. Standard theory shows that's all the solutions you can have--just two dimensions of them--and then the rest of the analysis goes through as shown here. – whuber May 15 '13 at 12:39
• @whuber thanks, yes, I was hoping OP will catch the general idea and derives details that he needs on his own. You are right - I should have mentioned that this is basically a perturbed form up to linear term - I'll put it in the text. – Vitaliy Kaurov May 15 '13 at 15:16

If we multiply through the second equation by t, NDSolve can find a solution using "EquationSimplification" -> "Residual":

sol = NDSolve[{
{x'[t] == y[t], t y'[t] == -y[t] - t Exp[x[t]]},
{x == 1, y == 0}},
{x, y}, {t, 0, 1}, Method -> {"EquationSimplification" -> "Residual"}]

Plot[{x[t], y[t]} /. sol // Flatten // Evaluate, {t, 0, 1}] 