2
$\begingroup$

I would like to plot the function v[r], and resolve this differential equation in spherical simmetry:

ClearAll[v, r, a, b, ϵ]

ode1 = {v'[r] == v[r]* ({2 a^2/r} - {b/{r^2}})/({{v[r]}^2 - a^2}), v[ϵ] == ϵ};

With[{a = 10^22, b = 10^6}, 
  Block[{ϵ = $MachineEpsilon}, 
    sol = NDSolve[ode1, v[r], {r, .1, 10^50}, SolveDelayed -> True]]];

Plot[Evaluate[v[r] \.sol], {r, .1, 10^50}, PlotRange -> All]

and the out is:

NDSolve::nlnum: The function value {{{0. -(2.22045*10^-16 (9.0072*10^15 Power[<<2>>]-2.02824*10^31 b))/(4.93038*10^-32+Times[<<2>>])}}} is not a list of numbers with dimensions {1} at {r,v[r],(v^′)[r]} = {2.22045*10^-16,2.22045*10^-16,0.}

and

NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions

I've looked for the errors and then I've tried to rewrite the equation, but the output is still like above.

$\endgroup$
  • 1
    $\begingroup$ Do not use { and } in place of parentheses. They have different meanings. $\endgroup$ – bbgodfrey Oct 24 '15 at 16:38
  • $\begingroup$ @bbgodfrey Yes,you're right but the errors remain the same $\endgroup$ – Lbk Oct 24 '15 at 16:55
  • $\begingroup$ b is undefined...Use Block instead of With, maybe. (You can see the symbol b in the error message.) -- Consider what happens here: foo = b; With[{b = 1}, foo] $\endgroup$ – Michael E2 Oct 24 '15 at 16:58
  • $\begingroup$ Also, you can use DSolve instead of NDSolve. $\endgroup$ – bbgodfrey Oct 24 '15 at 17:06
1
$\begingroup$

Starting with your expression let's change the inner curly brackets

ode1 = {v'[r] == v[r]*({2 a^2/r} - {b/{r^2}})/({{v[r]}^2 - a^2}), v[ϵ] == ϵ};

to parenthesize

ode1 = {v'[r] == v[r]*((2 a^2/r) - (b/(r^2)))/(((v[r])^2 - a^2)), v[ϵ] == ϵ}

Mathematica graphics

There are more parenthesis present than needed. For example, Power binds more tightly than Times. So it could be simplified a bit to:

ode1 = {v'[r] == v[r]*(2 a^2/r - b/r^2)/(v[r]^2 - a^2), v[ϵ] == ϵ}

I think it is possible to use Block to achieve the result but I'm simply going to create a new equation and substitute in the numerical values you selected.

ode2 = ode1 /. {a -> 10^22, b -> 10^6, ϵ -> $MachineEpsilon}

Mathematica graphics

Now execute NDsolve.

sol = NDSolve[ode2, v[r], {r, .1, 10^50}]

Mathematica graphics

You had a small mistake in your plot. You used \. when it should have been /.

I found the plot uninteresting over the entire range (I leave it to you) but near the origin:

Plot[v[r] /. sol, {r, 0.1, 10}, PlotRange -> All]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Note that the near-the-origin plot is still basically 0 throughout. Everything's of the order $10^{-31}$. $\endgroup$ – Patrick Stevens Oct 25 '15 at 15:47
  • $\begingroup$ @bbgodfrey : It is too slow to resolve the equation through DSolve ! @Jack LaVigne : I would to plot the same differential equation but with the following variables and axes: x= s= r/10^10, y=z(r)=v(r)/10 .Furthermore is there a kind of "time-shortcut" to resolve it through DSolve ? $\endgroup$ – Lbk Oct 26 '15 at 16:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.