Solve Nonlinear Stiff ODE in Mathematica (NDSolve)

Click For Summary

Discussion Overview

The discussion revolves around solving a nonlinear stiff ordinary differential equation (ODE) using Mathematica's NDSolve function. Participants explore issues related to numerical solutions, plotting results, and compatibility with other Mathematica functions like Manipulate.

Discussion Character

  • Technical explanation
  • Debate/contested
  • Homework-related

Main Points Raised

  • One participant attempts to solve the equation y''*y=-0.5*y'^2 using NDSolve but encounters an empty plot, questioning whether NDSolve can handle singular points.
  • Another participant suggests that the solution is complex and that the initial condition leads to complex values, recommending the use of the real root instead.
  • A participant inquires about the compatibility of NDSolve with Manipulate, expressing difficulty in combining the two.
  • One response emphasizes the importance of checking each component of the code and questions the appropriateness of the plotting range given the original solution's limits.
  • Another participant advises against using the % operator in Mathematica, citing potential confusion in later evaluations.

Areas of Agreement / Disagreement

Participants express differing views on the issues encountered with NDSolve and its syntax. There is no consensus on the best approach to resolve the problems, and multiple perspectives on how to handle the equations and plotting exist.

Contextual Notes

Limitations include potential misunderstandings of Mathematica's syntax and the complexity of NDSolve, which may lead to errors in implementation. The discussion also highlights the challenge of working with complex solutions and singular points.

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.
 

Similar threads

  • · Replies 24 ·
Replies
24
Views
4K
  • · Replies 1 ·
Replies
1
Views
925
  • · Replies 3 ·
Replies
3
Views
3K
  • · Replies 1 ·
Replies
1
Views
3K
  • · Replies 5 ·
Replies
5
Views
2K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 1 ·
Replies
1
Views
1K
  • · Replies 1 ·
Replies
1
Views
7K
  • · Replies 5 ·
Replies
5
Views
3K
  • · Replies 3 ·
Replies
3
Views
4K