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

Differential eq with square root numerical

  1. Aug 18, 2011 #1

    I have to solve a system of two coupled nonlinear ordinary differential equations. In the function that defines the system there are terms with square root of the dependant variables. I tried to solve the system numerically but in the nmerical solution I get violent oscillations. Then I realized that the same problem appears in the numerical solution of this very simple exactly solvable differential equation: x'(t)=(1-x(t))^(1/2) with initial condition x(2.1)=0.9975. The exact solution is x(t)=(1/4)(4t-t^2).

    It seems to me that the problem is the square root (1-x(t))^(1/2). When the numerical solution approach an extremum, then x'(t) and (1-x(t))^(1/2) both approach zero. But then (1-x(t)) also approaches zero, and in fact (1-x(t)) takes negative values, that makes the square root imaginary.

    I explore the situation with the fourth-order Runge-Kutta method with step h=0.2 and in the firt iteration I find the complex value:

    x(2.1)=1.00463+ 0.0042963 I

    The same happens with the Euler method and the MidPoint Method

    Any hint of how to get the correct numerical solution of this equation???

    Thanks a lot!
  2. jcsd
  3. Aug 18, 2011 #2


    User Avatar
    Homework Helper

    Can you write the system down.
  4. Aug 18, 2011 #3

    Initially I want to get the numerical solution of this very simple system.

    x'(t)=(1-x(t))^(1/2) with initial condition x(2.1)=0.9975.

    Also I note that if I start with initial conditions in the maximum x(2)=1, the RK4 gives contant solution x(t)=1. That is the reason I have to start with x(2.1)=0.9975.
  5. Aug 18, 2011 #4


    User Avatar
    Homework Helper

    You know that you can solve that analytically right? Another way to tackle this problem is by Newton's method, you define your points [itex]t_{i}[/itex], then you write down your equation at the point [itex]t_{i}[/itex], this then will lead to you computing the Jacobian and then it should only take a few iterations to get the answer.
  6. Aug 18, 2011 #5
    Yes! I know the exact solution. But I need to get the numerical solution of the differential equation. It is a warm up to a harder equation with similar problems.
  7. Aug 18, 2011 #6


    User Avatar
    Homework Helper

    I found that Newtons method works very well for this kind of thing but you have to have a reasonable first guess though, or at least something random that the iteration can work with.
  8. Aug 19, 2011 #7
    You could hard-code the algorithm to follow the correct path over the underlying multisheet of the square root function, or eliminate the multi-function from the differential equation:





    Now, in place of your original DE, use:

    [tex]x''=-1/2;\quad x(t_0)=x_0,\quad x'(t_0)=\sqrt{1-x}[/tex]
  9. Aug 19, 2011 #8


    User Avatar
    Homework Helper

    Jackmell, I think that this is a little specialised for the problem in hand (although it is a rather nice solution), it is just a warm up problem for him, I think he is asking for a more general method.
  10. Aug 19, 2011 #9
    Let's do the easy one first:

    The code below is the first approach in Mathematica using a modified Euler-method which maintains the correct path over the multisheet (Riemann surface of the complex square root).

    Code (Text):
    platch = True;
    xold = 0;
    thet = 0;
    a = 1;
    dh = 0.01;
    myvals = Table[temp = 1 - xold;
         If[platch == True && temp < 0,
          {xold = 2 - xold; platch = False; a = -1; }];
         xnew = xold + a*Sqrt[1 - xold]*dh;
         xold = xnew; {thet += dh, xnew}, {n, 1, 500}];
    ListPlot[myvals, PlotStyle -> Red, Joined -> True]
    And that gives the desired answer. Anyway, I think the problem with a standard numerical analysis of this DE, is that the algorithm is passing through the branch-point where the two branches of the root meet and since delta h is finite, it over-shoots the branch-point landing on a surface which has an imaginary component. In the algorithm above, I force the algorithm to follow a path which is always real. Maybe though I don't have this quite presented right but that's what I'm sticking with for now.

    How about posting the harder one xapiens.
    Last edited: Aug 19, 2011
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook