Message Boards Message Boards

Strange numerical problem with fairly simple symbolic equation.

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

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

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

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

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

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

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

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
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract