Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving non-linear ODE

  1. Mar 23, 2015 #1

    joshmccraney

    User Avatar
    Gold Member

    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)?
     
  2. jcsd
  3. Mar 23, 2015 #2
    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]
     
  4. Mar 23, 2015 #3

    joshmccraney

    User Avatar
    Gold Member

    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!
     
  5. Mar 23, 2015 #4
    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 whats going wrong near [itex] y=0[/itex]. You might have to adjust the domain to ensure that [itex] y [/itex] does not cross zero.
     
  6. Mar 23, 2015 #5

    joshmccraney

    User Avatar
    Gold Member

    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.
     
  7. Mar 23, 2015 #6
    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.
     
  8. Mar 24, 2015 #7

    joshmccraney

    User Avatar
    Gold Member

    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
     
  9. Mar 24, 2015 #8
    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.
     
  10. Mar 24, 2015 #9

    joshmccraney

    User Avatar
    Gold Member

    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?
     
    Last edited: Mar 24, 2015
  11. Mar 24, 2015 #10
    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}]
     
  12. Mar 24, 2015 #11

    joshmccraney

    User Avatar
    Gold Member

    Perfect! Thanks so much! Just curious, but why is it y[x]/.s? Why divide by .s?
     
  13. Mar 24, 2015 #12
    I think " /. " is the command to replace something with a predefined object, not exactly sure though.
     
  14. Mar 24, 2015 #13

    joshmccraney

    User Avatar
    Gold Member

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook