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$
4
  • 1
    $\begingroup$ Do not use { and } in place of parentheses. They have different meanings. $\endgroup$
    – bbgodfrey
    Commented Oct 24, 2015 at 16:38
  • $\begingroup$ @bbgodfrey Yes,you're right but the errors remain the same $\endgroup$
    – Lbk
    Commented Oct 24, 2015 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
    Commented Oct 24, 2015 at 16:58
  • $\begingroup$ Also, you can use DSolve instead of NDSolve. $\endgroup$
    – bbgodfrey
    Commented Oct 24, 2015 at 17:06

1 Answer 1

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$
2
  • $\begingroup$ Note that the near-the-origin plot is still basically 0 throughout. Everything's of the order $10^{-31}$. $\endgroup$ Commented Oct 25, 2015 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
    Commented Oct 26, 2015 at 16:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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