Solve Nonlinear Stiff ODE in Mathematica (NDSolve)

Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
5 replies · 6K views
Udi
Messages
3
Reaction score
0
Hello,

I'm trying to (numerically) solve the equation y''*y=-0.5*y'^2 in Mathematica.
I know there's an analytic solution (and I know how to calculate it), but I want to modify this equation and thus need to verify that the numerical solution for the original equation is exact.

I'm using NDSolve with the following syntax:

s = NDSolve[{y''[x] y[x] == -0.5 y'[x] y'[x],
y[0] == -75*((-4/675)^(1/3)), y'[0] == (-4/675)^(1/3)},
y, {x, 0, 48}]
Plot[(y /. Flatten[%])[x], {x, 0, 48}]

I know that this set of initial conditions have a well defined solution (solved analytically and numerically in Matlab) yet I get an empty plot in Mathematica. What is wrong? Could NDSolve handle equations with a singular point (here [50, 0])?
 
Physics news on Phys.org


The solution is complex, and Plot refuses to plot complex numbers. The problem is that -75*((-4/675)^(1/3)), when reduced to a number, is evaluated to -6.78604 - 11.7538 I. (When Plot gives you a blank, it's always a good idea to evaluate your function for a specific value or two.) You want the real root. I would just take the minus sign out of the cube root.

Code:
s2 = NDSolve[{y''[x] y[x] == -0.5 y'[x] y'[x], 
   y[0] == 75*((4/675)^(1/3)), y'[0] == -(4/675)^(1/3)}, 
  y, {x, 0, 48}]
Plot[y[x] /. s2, {x, 0, 48}]

This works fine.
 


Thanks!

By the way, is NDSolve compatible with Manipulate?
I've tried combining the two with little success:

Manipulate[
Plot[NDSolve[{y''[x]*(y[x] + eps/y[x]) == -0.5 y'[x] y'[x],
y[0] == 75*((4/675)^(1/3)), y'[0] == -(4/675)^(1/3)},
y, {x, 0, 48}], {x, 0, 100}], {eps, 0, 1}]

Again, this could be just a bad syntax problem.
 


Yes, NDSolve is compatible with Manipulate.

When something like this fails, you need to look at each of the pieces and see whether it's doing what you want. In this case I think you would know what the problem is if you just stared at it a little, because you used the right syntax earlier.

What, exactly, does NDSolve return? (A list of lists of replacement rules.) Can you plot that? (No.) Also, do you really think it's a good idea to plot the function from 0 to 100 when you only solved the DE from 0 to 48?

C'mon, man. Think!
 


I'm sorry, I'm new to mathematica, and don't know the syntax by heart yet (and what kind of variable each function returns).
Iv'e changed the boundaries to 0-100 in both case. Still doesn't work.
I know, I still need to change the input as to plot y instead of NDSolve. Not that I know how to this inside another function, but I'll think about it.
 


Yeah, it's tough, and NDSolve is not a good place to start. It's just about the most complex and finicky function in Mathematica.

It's best to build of complex expressions in pieces. Don't go for broke. Write some little part of it, evaluate it, and see what you get. Edit it if it's not what you expect. Then write another little part and evaluate that. Then combine the two and evaluate that. And so on.

One other thing I would advise: Never, never, never use %. I know Wolfram uses it in lots of examples, but it will confuse the Hell out of you later, because its value is the last thing that was evaluated, not the expression before the one you're now working on.