Re: NDSolve with square root

• To: mathgroup at smc.vnet.net
• Subject: [mg120999] Re: NDSolve with square root
• From: Peter Pein <petsie at dordos.net>
• Date: Sun, 21 Aug 2011 05:31:49 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• References: <j2o1n1\$4tf\$1@smc.vnet.net>

```Hi William,

I get two solutions, but there is just a linear shift along the t-axis:

In[1]:= Clear[x1, x2];
InputForm[{x1, x2} = x /. DSolve[{Derivative[1][x][t] == Sqrt[1 - x[t]],
x[21/10] == 399/400}, x, t]]

Out[2]//InputForm=
{Function[{t}, (-21 + 110*t - 25*t^2)/100], Function[{t}, (4*t - t^2)/4]}

In[3]:= Reduce[ForAll[t, x1[t] == x2[a + b*t]] && b >= 0, {a, b}]//InputForm

Out[3]//InputForm=
a == -1/5 && b == 1

check:
In[4]:= Expand[x2[x - 1/5] - x1[x]]
Out[4]= 0

If you insist on floats, try
{x1, x2} = Composition[N, x] /. DSolve[...]

hth,
Peter

Am 20.08.2011 12:19, schrieb william cuervo:
> 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] == Sqrt[1 - x[t]] , x[2.1] == 0.9975
>
> The exact solution is x(t)=(1/4)(4t-t^2).
>
> I tried:
> NDSolve[{x'[t] == Sqrt[1 - x[t]] , x[2.1] == 0.9975`}, x, {t, -5, 5},
> MaxSteps ->  4000000]
>
> using ExplicitEuler, ExplicitRungeKutta and MidPoint Methods.
>
> 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!

```

• Prev by Date: Re: Display BarChart ticks in AccountingForm?
• Next by Date: Re: More robust pattern based replacement rules?
• Previous by thread: Re: NDSolve with square root
• Next by thread: Getting 3d printable models out of 3d graphics in Mathematica