Message Boards Message Boards

0
|
5587 Views
|
10 Replies
|
1 Total Likes
View groups...
Share
Share this post:

Solve second order stiff ODE with boundary conditions?

Posted 8 years ago

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

POSTED BY: satoru none
10 Replies

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!

plot of solution

POSTED BY: Kay Herbert

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

POSTED BY: Hans Dolhaine

Trying to solve second order stiff ode with boundary conditions

POSTED BY: Kay Herbert

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]]}}
POSTED BY: Kay Herbert

@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 ].

POSTED BY: Hans Dolhaine
Posted 8 years ago

Any ideas what to do? Any one know good external packages to solve stiff systems?

POSTED BY: satoru none

@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
POSTED BY: Hans Dolhaine

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

POSTED BY: Hans Dolhaine

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?

POSTED BY: Hans Dolhaine
Posted 8 years ago

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.

POSTED BY: satoru none
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