How do I numerically solve a non-linear ODE in Mathematica?

[member 428835]

Hi PF!

I am wondering if any of you have experience numerically solving second order ODE's? Basically, I'm trying to solve one and am trying to do it numerically in mathematica. Can anyone help? For those curious, the equation is ##y y'' + 2y'^2 +xy' = 0## where ##y## is a function of ##x##.

Thanks so much!

PS is anyone else unable to see latex on physics forums (mine's simply showing all statements literally)?

 

[HomogenousCow]

I would just use NDSolve.
The code is
NDSolve[{y[x] y''[x]+2(y'[x])^2+x(y'[x])==0,y'[x1]==a,y[x1]==b},y[x],{x,x1,x2}]

I've assumed you want both initial conditions at the point x1, and the equation solved on the interval [x1,x2]

 

[member 428835]

When I try to put that in I receive an error message, something about a singularity and 1/0. My input is below

NDSolve[{x Derivative[1][y][x] + 2 Derivative[1][y][x]^2 + y[x] (y^\[Prime]\[Prime])[x] == 0, Derivative[1][y][1] == -(1/2), y[1] == 0}, y[x], {x, 0, 1}]

the boundary conditions are ##y(1)=0## and ##y'(1)=1/2##.

Please let me know how you would proceed, and thanks so much for helping me!

 

[the_wolfman]

Your differential equation is singular where [itex]y=0[/itex]. Here the differential equation is essentially [itex] \left(0\right)y'' +2y'^2+xy'=0[/itex]. You'll note that here there is no constraint on [itex] y'' [/itex] because it is multiplied by zero. Here normal existence and uniqueness theorems fail. I suspect that your initial conditions are not uniquely specifying the solution. They are redundant. Try solving the equation with initial conditions [itex] y \ne 0 [/itex]. This might help show what's going wrong near [itex] y=0[/itex]. You might have to adjust the domain to ensure that [itex] y [/itex] does not cross zero.

 

[member 428835]

Thanks for taking the time to reply. So I used the only other boundary condition I know of (although this one evidently was numerically found) as ##y(0)=.345...## and I'm still getting an error. But from the text it sounds like ##y(0)=.345...## was found only after solving the initial ODE with the above listed boundary conditions.

Any other idea of how to solve, or can mathematica not do this? Again, I'd like to try and numerically solve using the first two boundary conditions.

 

[the_wolfman]

I'm not sure I understand your last comment. I encourage you to try a bunch of different initial y and y', use a different domain, and a different initial x. Just pick random numbers. Avoid y=0.

There are two things going on. First numerically, I suspect that Mathematica might be running into errors associated with dividing by zero when y=0. I'm not too familiar with Mathematica, but I suspect its trying to reformulate the equation into something akin to [itex] y'' =\frac{1}{y} \times f(x,y,y')[/itex]. You could try modifying the equation to [itex](y+\epsilon)y''+2y'^2+xy'=0 [/itex] where [itex] \epsilon [/itex] is a small positive number. See what happens.

I also suspect that there is a fundamental error in your initial conditions. I suspect that there are multiple solutions that pass through the point [itex] y(1)=0, y'(1)=1/2 [/itex]. Singular differential equations often have this problem. If this is true, then your problem as stated is not well posed. There will be many solutions to the differential equation with the given initial conditions. You may be able to trick Mathematica into giving you one solution. But you really need to understand the behavior near the singularity before you trust any solution. This is why I encourage you to try many different initial conditions.

 

[member 428835]

OK great, I've tried setting the zero boundary condition to a number close to zero and am getting great results, so thanks!

However, can you explain what you meant when you said that

normal existence and uniqueness theorems fail.

 

[HomogenousCow]

OK great, I've tried setting the zero boundary condition to a number close to zero and am getting great results, so thanks!

However, can you explain what you meant when you said that

He means that when a DE is singular at a point, boundary conditions applied there may not lead to a unique solution, or any solution.

 

[member 428835]

Ok, thanks for explaining that. So since, as was previously stated, mathematica is probably isolating ##y''## then because ##y## is in the denominator we have a problem when ##y=0## right?

Also, I was hoping I could integrate my numerical solution I received in mathematica. Specifically, I am trying to calculate ##\int_0^1 y^2 \, dx##. Do you know how to do this? I get a result if I type
Integrate[(y[x] /. t), {x, 0, 1}]
where I used ndsolve over the above ode, but for some reason if i type
Integrate[(y[x] /. t)^2, {x, 0, 1}]
i get no answer. Any ideas?

 

[HomogenousCow]

Ok, thanks for explaining that. So since, as was previously stated, mathematica is probably isolating ##y''## then because ##y## is in the denominator we have a problem when ##y=0## right?

Also, I was hoping I could integrate my numerical solution I received in mathematica. Do you know how to do this? I can give the technique used to solve it if that helps?

Sure that's easy to do.
First call your solution something, to do this you want to give your NDSolve[..] a name.
For example,
s=NDSolve[...] (where ... represents all the stuff you have inside the function)
Then after you shift-enter, you can use the function in any other operation by typing "y[x]/.s" in place of the function
If you want to integrate this function from {a,b} (And make sure the function has been solved on this interval), use
NIntegrate[y[x]/.s, {x,a,b}]

 

[member 428835]

Perfect! Thanks so much! Just curious, but why is it y[x]/.s? Why divide by .s?

 

[HomogenousCow]

Perfect! Thanks so much! Just curious, but why is it y[x]/.s? Why divide by .s?

I think " /. " is the command to replace something with a predefined object, not exactly sure though.

 

[member 428835]

thanks!