Differential eq with square root numerical

[xapiens]

Hi!

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!

 

[hunt_mat]

Can you write the system down.

 

[xapiens]

Hi!

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.

 

[hunt_mat]

You know that you can solve that analytically right? Another way to tackle this problem is by Newton's method, you define your points t_{i}, then you write down your equation at the point t_{i}, this then will lead to you computing the Jacobian and then it should only take a few iterations to get the answer.

 

[xapiens]

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.

 

[hunt_mat]

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.

 

[jackmell]

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

Thanks a lot!

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:

x'=\sqrt{1-x}

x''=1/2(1-x)^{-1/2}(-x')

x''=-1/2\frac{x'}{\sqrt{1-x}}

x''=-1/2

Now, in place of your original DE, use:

x''=-1/2;\quad x(t_0)=x_0,\quad x'(t_0)=\sqrt{1-x}

 

[hunt_mat]

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.

 

[jackmell]

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.

Let's do the easy one first:

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.

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).

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.