I am trying to handle numericaly well known exactly solvable problem:
\[Lambda] = 0.7
v = 100
NDSolve[{f''[x] == \[Lambda]^2/2 f[x]*((f[x])^2 - v^2), f[-1000] == v,
f[1000] == -v}, f[x], {x, -1000, 1000}]
I use 1000 instead of infinity.
This must be kink solution, but for all methods it fails with NDSolve::ndsz: At x == -999.964, step size is effectively zero; singularity or stiff system suspected. >>
Sometimes rescaling of the problem will do a good job. I didn’t look too much in your system, but with my version (Mma 7) and a modified equation things seem to work:
\[Lambda] = 0.7/1000
v = 100
fff = f[x] /. NDSolve[{f''[x] == \[Lambda]^2/2 f[x]*((100 f[x])^2 - v^2), f[-1] == v, f[1] == -v}, f[x], {x, -1, 1}][[1, 1]]
Plot[fff, {x, -1, 1}]
To be discussed: is that rescaling (which in effect is a change of coordinates) valid, resp. does it give an equivalent system?
Thx for reply! Rescaling is a good idea, but i think u did it wrong. Your solution dont contains horizontal plateau (like Tanh[x] - exact solution of this system) near the borders of integration - reason of equation stiffness. Whats why u was able to get it. My try is \[Lambda] = 0.7 v = 100 \[Mu] = 10 fff = NDSolve[{f''[ x] == \[Mu]^2 (\[Lambda]^2/2) f[x]*((f[x])^2 - v^2), f[-1] == v, f[1] == -v}, f\[x\], {x, -1, 1}, MaxSteps -> 1000000]
Here mu is rescaling factor, tried from 1000 to 1 and allways failed with stiffness report.
Hello,
but i think u did it wrong
of course you are right. I did my rescaling “quick and dirty”, just to give an idea. I know a couple of examples where NDSolve “crashed” and a rescaling just made it. I tried it once more:
fT[x_] := A \[Phi][ B x]
\[Lambda] = .7; v = 100;
eq = f''[x] == \[Lambda]^2/2 f[x] (f[x]^2 - v^2)
In[9]:= eq2 = #/(A B^2) & /@ eq1 /. x -> z/B // Expand
Out[9]=
\!\(\*SuperscriptBox["\[Phi]", "\[Prime]\[Prime]",MultilineFunction->None]\)[z] == -((2450. \[Phi][z])/B^2) + (
0.245 A^2 \[Phi][z]^3)/B^2
and played around with A and B, but the stiffness error persisted.
Then I broke down your diff.eq. of second order into two of first order (but I am quite sure that is exactly what NDSolve does internally), but everything stayed as before. Seems you have a system which can’t be attacked by NDSolve…
Any ideas what to do? Any one know good external packages to solve stiff systems?
Potentially you can find an analytical solution to your equation.
Integrating your equation you can find:
f’[x]^2 == lambda^2 (1/2 f^4 - v^2 f^2 + 1/2 v^4)
or
f’[x] == - lambda Sqrt[1/2 f^4 -v^2 f^2 + 1/2 v^4]
You can now numerically integrate or solve analytically the first order equation (which is messy)
In[2]:= DSolve[
f'[x] == - lambda Sqrt[1/2 f[x]^4 - v^2 f[x]^2 + 1/2 v^4], f[x], x]
Out[2]= {{f[x] ->
InverseFunction[-((ArcTanh[#1/v] (-v^2 + #1^2))/(
v Sqrt[(-v^2 + #1^2)^2])) &][-((lambda x)/Sqrt[2]) + C[1]]}}
@Kay
That is a very interesting idea.
But I get a different integrand:
We had
eq1 = D[f[x], {x, 2}] == \[Lambda]^2/2 f[x] (f[x]^2 - v^2) // Expand
which by integration gives
eq2 = 1/2 D[f[x], x]^2 == \[Lambda]^2/2 (-(1/2) v^2 f[x]^2 + f[x]^4/4)
Let’s have a look:
D[#, x] & /@ eq2 // Simplify
Collect[Solve[%, f''[x]], \[Lambda]]
which gives indeed an equivalent to eq1.
But
Solve[eq2, D[f[x], x]]
seems to give a different expression for f ’ [ x ] than given in your post.
And finally: unfortunately we have no information about f ’ [ -1000 ].
@satoru
Well, things really seem to be quite complicated. I think it is astonishing that NDSolve crashes all the time. So I thought what one could do. Just giving your diff.equation leave no options. So I decided to treat it somewhat more flexible. Rewriting your diff.eq, gives a system of to equations of 1st order:
\[Lambda]n = .7; vn = 100;
sys1 = {
f'[x] == u[x],
u'[x] == \[Lambda]n f[x] (f[x]^2 - vn^2),
f[-1000.] == vn,
u[-1000.] == A};
The idea is, if this system could be integrated, to vary the value A and look for the value of f [ 1000 ], hopefully finding an A that “hits” the right value.
So I tried (note the huge value for MaxSteps)
Clear[sol1]
sol1[An_] := Module[{},
sys = sys1 /. A -> An;
(*Print[sys];*)
NDSolve[sys, {f[x], u[x]}, {x, -1000., 1010.}, MaxSteps -> 10^7] //
Flatten]
and
In[122]:= ff = f[x] /. sol1[-5. 10^-6];
ff /. x -> 1000
Out[123]= -97.4093
That is almost the -100 you wanted to have. The result is a fast varying oscillation (note that i plotted it only i a small window of width ten):
Plot[ff, {x, -1000, -990}]
Is such an oscillation what you expected in your problem?
BUT : the behavior of that system seems to be quite insensitive with respect to the value of A (if it does not crash, which occurs as well). Note that I enhanced A by a factor of 1000
In[129]:= sol1[-.5 10^-3];
ff /. x -> 1000
Out[130]= -97.4093
Trying to solve second order stiff ode with boundary conditions
Hi,
Your f goes to fixed values as x goes to +/- infinity consequently f’ needs to go to 0 as x goes to +/- infinity.
I did the integration in my head so I think I was off by a factor of 2. There is an integration constant which I set so that f’ vanishes at +/- infinity:
In[14]:= eq2 =
D[f[x], x]^2 == \[Lambda]^2 (-1/2 v^2 f[x]^2 + 1/4 f[x]^4 + 1/4 v^4)
Out[14]= Derivative[1][f][
x]^2 == \[Lambda]^2 (v^4/4 - 1/2 v^2 f[x]^2 + f[x]^4/4)
In[15]:= D[#, x] & /@ eq2 //
Simplify Collect[Solve[%, f''[x]], \[Lambda]]
Out[15]= {}[
2 Derivative[1][f][x] (f^\[Prime]\[Prime])[
x] == \[Lambda]^2 (-v^2 f[x] Derivative[1][f][x] +
f[x]^3 Derivative[1][f][x])]
Which gives your equation 1
I then take the square root of each side picking the correct branch
In[74]:= eq3 =
f'[x] == -\[Lambda] Sqrt[-1/2 v^2 f[x]^2 + 1/4 f[x]^4 + 1/4 v^4]
Out[74]= Derivative[1][f][x] == -\[Lambda] Sqrt[
v^4/4 - 1/2 v^2 f[x]^2 + f[x]^4/4]
and then I integrated to find x:
In[268]:= Integrate[-1/(\[Lambda] Sqrt[-1/2 v^2 f^2 + 1/4 f^4 +
1/4 v^4]), f]
Out[268]= (2 (f^2 - v^2) ArcTanh[f/
v])/(v Sqrt[(f^2 - v^2)^2] \[Lambda])
now simplify and you find:
In[273]:= Clear[v, \[Lambda]]
In[274]:= f[x_] := -v Tanh[\[Lambda] v x /2]
ok now I check if I did my algebra right this time:
In[275]:= eq1
Out[275]=
1/2 v^3 \[Lambda]^2 Sech[(v x \[Lambda])/2]^2 Tanh[(v x \[Lambda])/
2] == 1/2 v^3 \[Lambda]^2 Tanh[(v x \[Lambda])/2] -
1/2 v^3 \[Lambda]^2 Tanh[(v x \[Lambda])/2]^3
In[276]:= Simplify[%]
Out[276]= True
eureka!
@Kay
That is really cool. I overlooked the integration-constant. Without it the integration yields strange looking complex entities.
The remaining question is: what happens with NDSolve?
