Message Boards Message Boards

Strange numerical problem with fairly simple symbolic equation.

xort dsc
Posted 10 years ago

Hi,

I have a strange problem with a fairly simple function. I define the following:

In[1]:= u = Tanh[x/Sqrt[2]]
Out[1]= Tanh[x/Sqrt[2]]

In[2]:= f = Laplacian[u, {x}] + (1 - u^2)* u
Out[2]= -Sech[x/Sqrt[2]]^2 Tanh[x/Sqrt[2]] + Tanh[x/Sqrt[2]] (1 - Tanh[x/Sqrt[2]]^2)

Looks not too complicated, does it ? Now to the weired thing: I believe that u=Tanh[x/Sqrt[2]] is a solution to Laplacian[u, {x}] + (1 - u^2)* u == 0 using the boundary condition f[-inf]==-1 && f[inf]==1. Mathematica doesn't seem to be able to deliver me any solution, though. Well, I might be wrong with my guessed solution, so I just plugged it into the equation to see how far off I am.

Looking at the plot of f(x) from x=-4 to x=4 reveals something strange, though:

Weired plot of f from x=-4 to x=4

So that plot looks way too noisy and within tiny numbers that I suspect a numerical problem here (and tanh(x/sqrt(2)) actually being a solution). Is that possible ? Why is that ? Could that also be the problem why Mathematica is unable to solve for u with the said boundary condition ?

Also I tried

In[5]:= Simplify[f]
Out[5]= 0

so Mathematica DOES understand that it is zero. Why can't I then solve for u ? How would I do that properly ?
POSTED BY: xort dsc
15 Replies
Sort By:

For the boundary condition one gets an error:

DSolve[{D[u[x], {x, 2}] + (1 - u[x]^2) u[x] == 0, u[-Infinity] == -1, u[Infinity] == 1}, u, x]
Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>
DSolve::bvlim: For some branches of the general solution, unable to compute the limit at the given points. Some of the solutions may be lost. >>
Solve::incnst: Inconsistent or redundant transcendental equation. After reduction, the bad equation is 1-2 C[1] == 0. >>

the equation has a solution:

In[21]:= DSolve[D[u[x], {x, 2}] + (1 - u[x]^2) u[x] == 0, u, x]

During evaluation of In[21]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

Out[21]= {{u -> 
   Function[{x}, 
    I Sqrt[-(1/(1 - Sqrt[1 - 2 C[1]]))]
       JacobiSN[Sqrt[
       x^2 + x^2 Sqrt[1 - 2 C[1]] + 2 x C[2] + 
        2 x Sqrt[1 - 2 C[1]] C[2] + C[2]^2 + Sqrt[1 - 2 C[1]] C[2]^2]/
       Sqrt[2], (1 - Sqrt[1 - 2 C[1]])/(1 + Sqrt[1 - 2 C[1]])] - 
     I Sqrt[-(1/(1 - Sqrt[1 - 2 C[1]]))] Sqrt[1 - 2 C[1]]
       JacobiSN[Sqrt[
       x^2 + x^2 Sqrt[1 - 2 C[1]] + 2 x C[2] + 
        2 x Sqrt[1 - 2 C[1]] C[2] + C[2]^2 + Sqrt[1 - 2 C[1]] C[2]^2]/
       Sqrt[2], (1 - Sqrt[1 - 2 C[1]])/(
       1 + Sqrt[1 - 2 C[1]])]]}, {u -> 
   Function[{x}, -I Sqrt[-(1/(1 - Sqrt[1 - 2 C[1]]))]
       JacobiSN[Sqrt[
       x^2 + x^2 Sqrt[1 - 2 C[1]] + 2 x C[2] + 
        2 x Sqrt[1 - 2 C[1]] C[2] + C[2]^2 + Sqrt[1 - 2 C[1]] C[2]^2]/
       Sqrt[2], (1 - Sqrt[1 - 2 C[1]])/(1 + Sqrt[1 - 2 C[1]])] + 
     I Sqrt[-(1/(1 - Sqrt[1 - 2 C[1]]))] Sqrt[1 - 2 C[1]]
       JacobiSN[Sqrt[
       x^2 + x^2 Sqrt[1 - 2 C[1]] + 2 x C[2] + 
        2 x Sqrt[1 - 2 C[1]] C[2] + C[2]^2 + Sqrt[1 - 2 C[1]] C[2]^2]/
       Sqrt[2], (1 - Sqrt[1 - 2 C[1]])/(1 + Sqrt[1 - 2 C[1]])]]}}

The printing stuff below and around the $MachinePrecision means that it is not meaningful to print because it might or might not be zero, for Chop it's clearly zero.
POSTED BY: Udo Krause

Plot of necessity must evaluate numerically. That means there can be numerical error. The fact that you get values around $MachineEpsilon indicates this is in fact the case.

Sometimes one can fool Plot by adding a largish constant, say 10^4, in the computation and later subtracting it.
POSTED BY: Daniel Lichtblau 
u[x_] := Tanh[x/Sqrt[2]]
In[8]:= D[u[x], {x, 2}] + (1 - u[x]^2)*u[x]
Out[8]= -Sech[x/Sqrt[2]]^2 Tanh[x/Sqrt[2]] + 
 Tanh[x/Sqrt[2]] (1 - Tanh[x/Sqrt[2]]^2)  
In[9]:= Simplify[%] 
Out[9]= 0

Noisy graph just shows computational roundoff errors
POSTED BY: S M Blinder

There is a way

In[40]:= First[
          FullSimplify[
           ExpToTrig[u[y] /. DSolve[{D[u[x], {x, 2}] + (1 - u[x]^2) u[x] == 0, u[0] == 0, u'[0] == 1/Sqrt[2]}, u, x]]
          ]
         ]
Out[40]= Tanh[y/Sqrt[2]]

using conditions not so far away

In[41]:= DSolve[{D[u[x], {x, 2}] + (1 - u[x]^2) u[x] == 0, u[0] == 0, u'[0] == 1/Sqrt[2]}, u, x]
Out[41]= {{u -> Function[{x}, (-1 + E^(Sqrt[2] x))/(1 + E^(Sqrt[2] x))]}}
POSTED BY: Udo Krause

and the boundary condition at x^2+y^2+z^2==inf being {x,y,z}/Sqrt[x^2+y^2+z^2])

Why don't you adopt a radial symmetric condition? By the way possibly the above condition is for a scalar function $u(x,y,z)$ contradictory e.g. for {x, y, z} = {0, -Infinity, Infinity}. Not sure how this works for a vector valued function $\vec{u} = u_1(x,y,z) \vec{e_1} + u_2(x,y,z) \vec{e_2} + u_3(x,y,z) \vec{e_3}$.
POSTED BY: Udo Krause

In Mathematica 10 no cheating is needed, looking at the general solution one sees the C[1] and C[2] bringing the most significant simplifications:

In[44]:= Simplify[
    Quiet[(u[y] /. DSolve[D[u[x], {x, 2}] + (1 - u[x]^2) u[x] == 0, u, x])] /. {C[1] -> 1/2, C[2] -> 0}, 
    Assumptions -> y \[Element] Reals]
Out[44]= {-Tanh[Abs[y]/Sqrt[2]], Tanh[Abs[y]/Sqrt[2]]}

because Tanh is an odd function

In[45]:= Tanh[x] == -Tanh[-x]
Out[45]= True

this is your result. Seemingly the general algorithm has a problem to check that C[1] -> 1/2 und C[2] -> 0 will do the boundary conditions, saying

Solve::incnst: Inconsistent or redundant transcendental equation. After reduction, the bad equation is 1-2 C[1] == 0. >>

but this is in fact the good equation reducing the JacobiSN.
POSTED BY: Udo Krause
xort dsc
xort dsc
Posted 10 years ago

Thanks the replies. So is there no way Mathematica can solve this equation for the boundary condition and give me the tanh(x/sqrt(2)) ? Do I have to guess all my solutions for specific boundary problems ? ;)
POSTED BY: xort dsc
xort dsc
xort dsc
Posted 10 years ago

Even though I have luckily guessed the solution for the 1-dimensional case I was not so lucky with the 2- and 3-dimensional case (where u^2 becomes u.u and the boundary condition at x^2+y^2+z^2==inf being {x,y,z}/Sqrt[x^2+y^2+z^2]). Is there a way to solve this problem with Mathematica ?
POSTED BY: xort dsc
xort dsc
xort dsc
Posted 10 years ago

hmm, it doesn't work for me for some reason

In[5]:= First[FullSimplify[ExpToTrig[u[y] /. DSolve[{D[u[x], {x, 2}] + (1 - u[x]^2) u[x] == 0, u[0] == 0, u'[0] == 1/Sqrt[2]}, u, x]]]]
During evaluation of In[5]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>
During evaluation of In[5]:= DSolve::bvfail: For some branches of the general solution, unable to solve the conditions. >>
During evaluation of In[5]:= DSolve::bvfail: For some branches of the general solution, unable to solve the conditions. >>
Out[5]= y

In[6]:= DSolve[{D[u[x], {x, 2}] + (1 - u[x]^2) u[x] == 0, u[0] == 0, u'[0] == 1/Sqrt[2]}, u, x]
During evaluation of In[6]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>
During evaluation of In[6]:= DSolve::bvfail: For some branches of the general solution, unable to solve the conditions. >>
During evaluation of In[6]:= DSolve::bvfail: For some branches of the general solution, unable to solve the conditions. >>
Out[6]= {}

Any idea why that is ?

Also I don't understand how you got the new boundary condition from the original one.
POSTED BY: xort dsc
xort dsc
xort dsc
Posted 10 years ago

Yes, using spherical coordinates would certainly make sense as the solution must be spherically symmetric. I just wasn't sure how to do it with Mathematica (I'm still a beginner). And to reformulate the boundary condition: Basically the condition should be that vectors at r=inf should have a length=1 and point away from the origin. At the origin (r=0) the vector length=0. How could I formulate it for mathematica ? Though I'm still struggling with the 1D case (see comment to your next post).
POSTED BY: xort dsc

Any idea why that is ?

The Mathematica version is a candidate; here is Mathematica 10.0.2.0 Windows 7 64 Bit Home Premion in use.

Also I don't understand how you got the new boundary condition from the original one.

That's incredibly simple. You stated

In[2]:= Simplify[D[#, {x, 2}] + (1 - #^2) #] &[Tanh[x/Sqrt[2]]]
Out[2]= 0

So one uses the solution to compute other conditions

In[3]:= Limit[Tanh[x/Sqrt[2]], x -> 0]
Out[3]= 0

In[4]:= Limit[D[Tanh[x/Sqrt[2]], x] /. x -> y, y -> 0]
Out[4]= 1/Sqrt[2]

allowing to reproduce it

So is there no way Mathematica can solve this equation for the boundary condition and give me the tanh(x/sqrt(2)) ?

In[5]:= DSolve[{D[u[x], {x, 2}] + (1 - u[x]^2) u[x] == 0, u[0] == 0, u'[0] == 1/Sqrt[2]}, u, x]
Out[5]= {{u -> Function[{x}, (-1 + E^(Sqrt[2] x))/(1 + E^(Sqrt[2] x))]}}
POSTED BY: Udo Krause

with the 2- and 3-dimensional case (where u^2 becomes u.u ...)

In 2- and 3-D the equation $\begin{equation}\frac{d^2u(x)}{dx^2}+(1-u(x)^2) u(x)=0\end{equation}$ no obvious generalization for a vector valued function $\vec{u}(\vec{x})$ exists because the Laplace operator $\Delta=\nabla\cdot\nabla$ as a gradient of the divergence of a scalar function $u(\vec{x})$ will not happen to be applicable. If one applies the divergence to a vector valued function, a matrix valued function appears ...
POSTED BY: Udo Krause
xort dsc
xort dsc
Posted 10 years ago

Ah okay, but that's somewhat cheating as I have to know the solution in advance already. ;) If I wouldn't have had that lucky guess of tanh(x/sqrt(2)) in the first place there would be no way to derive the solution, i guess ?

(and i'm still on mathematica 9)
POSTED BY: xort dsc
xort dsc
xort dsc
Posted 10 years ago

hmm, i don't really understand what you are saying. But this equation is the Ginzburg-Landau equation which works in any dimension (so u yields a 3-vector for every given x, also a 3-vector). So for 3 dimensions it works on a 3d vector field. the laplacian of a 3d vector field is again a 3d vector field (and the second term is also a 3d vector field to be added). and i know from numerical experiments that there is a solution to it. so what's the problem here ?
POSTED BY: xort dsc

I was remembering only the scalar Laplacian but there is also the Vector Laplacian, you're right.
POSTED BY: Udo Krause
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract