Re: Difficulty in obtaining a numerical solution for a partial
- To: mathgroup at smc.vnet.net
- Subject: [mg114318] Re: Difficulty in obtaining a numerical solution for a partial
- From: Roland Franzius <roland.franzius at uos.de>
- Date: Wed, 1 Dec 2010 02:09:38 -0500 (EST)
- References: <id2eho$d2b$1@smc.vnet.net>
Am 30.11.2010 10:01, schrieb Dilip Ahalpara:
> Dear all,
>
> I have a partial differential equation (PDE), which when I try to solve
> using NDSolve (to get U(x,t) in a given range, say x in [0,20] and t in
> [0,10]),
> it blows up.
> The PDE for U(x,t) is,
> Ut + (2/U) Ux Uxx - Uxxx = 0
> where Ux=delU/delx, Ut=delU/delt, Uxx=del^2U/delx^2 etc.
>
> As initial condition, a profile U(x,0) is fed, which is obtained
> from a known function satisfying the PDE, obtained
> from
> (c0/2) * (Sech[Sqrt[c0]/2 * (x-x0-c0*t)])^2
> and the initial profile is placed around x0, with a height specified by c0.
> The boundary condition is set to match at x0Left and x0Right.
>
> To show that my Mathematica program has a correct logic,
> I first solve a PDE
> Ut + 6 U Ux + Uxxx = 0
> using the same initial and boundary conditions -- and it works
> nicely as expected.
>
> I include below the Mathematica program.
>
> ----------
> (* Define the differntial equations*)
> diffEqnLHS1 = D[u1[x, t], t] + 6*u1[x, t]*D[u1[x, t], x] + D[u1[x, t], {x,
> 3}]
> diffEqnLHS2 = D[u2[x, t], t] + D[u2[x, t], x] D[u2[x, t], {x, 2}]*2/u2[x,
> t] - D[u2[x, t], {x, 3}]
> template[x0_, c0_, x_, t_] := 0.5*c0*(Sech[0.5*Sqrt[c0]*(x - x0 - c0*t)])^2;
>
> cValue = 1.2;
> {x0Left, x0Right} = {0.0, 20.0}; (*Solution space*)
> x0Template = 5.0; (*Placing initial profile at x0*)
> c0Template = 3;(*with height c0*)
> initCondition = template[x0Template, c0Template, x, 0] +
> template[x0Template + (x0Right - x0Left), c0Template, x, 0] +
> template[x0Template - (x0Right - x0Left), c0Template, x, 0];
> Print[N[initCondition /. x -> x0Left, 10], " ", N[initCondition /.x ->
> x0Right, 10],
> " ", (initCondition /. x -> x0Left) == (initCondition /.x -> x0Right)];
>
> soln1 = NDSolve[diffEqnLHS1 == 0&& u1[x, 0] == initCondition&& u1[x0Left,
> t] == u1[x0Right, t],
> u1, {t, 0, 10}, {x, x0Left, x0Right}]; // Timing
> Animate[Plot[Evaluate[u1[x, t] /. soln1[[1]]], {x, x0Left, x0Right},
> PlotRange -> {0, 2}, GridLines -> Automatic], {t, 0, 10}]
>
> soln2 = NDSolve[diffEqnLHS2 == 0&& u2[x, 0] == initCondition&& u2[x0Left,
> t] == u2[x0Right, t],
> u2, {t, 0, 1}, {x, x0Left, x0Right}]; // Timing
> Animate[Plot[Evaluate[u2[x, t] /. soln2[[1]]], {x, x0Left, x0Right},
> PlotRange -> {0, 2}, GridLines -> Automatic], {t, 0, 1}]
> ----------
>
> As can be seen by running this code, the Animate for soln1 works
> quite ok -- it shows the initial profile propagating to the right
> (increasing x) without changing the shape of the profile.
> However when I try to solve another PDE of my interest,
> i.e. corresponding to soln2, NDSolve does not generate a
> well behaved solution.
>
> My questions:
> [1] What is going wrong while solving the second PDE?
Even if you avoid a zero in the second term 2/U by replacing it with
2/(1+0.01 U^2) the numerical solution becomes instable at t=0.01 and
quite suddenly obtains theets on straight sec^2 edges at a time the wave
shape has'nt moved at all in space. The problem is probably that at U_xx
=0 the term U_xxx generates exponential growth.
> [2] How to obtain a numerical solution for the second
> PDE correctly?
Generally speaking, the soliton solutions of nonlinear wave equations
represent an algebraically knotted system of functions, where the
solution can be exponentially small for a large finite x-t-area so that
numerically it seems not to be different in this area from a zero solution.
These areas then represent permanent sources of digital noise. That in
turn taken as a starting value generates exponential big lumps or kinks
somewhere far outside the area in time and space.
Using periodic boundary conditions in space you fold all the external
noise soliton maxima by the mirror principle into your x-intervall after
a short time.
Thats the reason why one uses hamiltonian methods to control numerically
as many constants of motion during the time evolution step evaluation as
possible, very similar like the methods using momentum, energy and
angular momentum control of subsets of particles in many particle
physics.
--
Roland Franzius