# Differential eq with square root numerical

1. Aug 18, 2011

### 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!

2. Aug 18, 2011

### hunt_mat

Can you write the system down.

3. Aug 18, 2011

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

4. Aug 18, 2011

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

5. Aug 18, 2011

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

6. Aug 18, 2011

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

7. Aug 19, 2011

### jackmell

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}$$

8. Aug 19, 2011

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

9. Aug 19, 2011

### jackmell

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