MathGroup Archive 2002

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Method of line

  • To: mathgroup at smc.vnet.net
  • Subject: [mg36215] Re: Method of line
  • From: Robert Knapp <rknapp at wolfram.com>
  • Date: Tue, 27 Aug 2002 02:07:22 -0400 (EDT)
  • Organization: Wolfram Research, Inc.
  • References: <akcot4$52h$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

amar ali wrote:
> Dear Mathgroup
> I have convert pde to ode like this
> 
> du_i/dt  =  1-4 u_i + .02 (u_i-1 ? 2 u_i + u i+1)/(delta (x))^3  + (u_i)^3 
> v_i
> 
> dv_i /dt  =  3 u_i + .02 (v_i-1 ? 2 v_i +  u i+1)/(delta (x))^3  - (u_I)^3 
> v_i
> 
> delta(x)=(i)/(N+1)
> 
> x_i=  (i)/(N+1)
> 
> Boundary condition
> u_0  (t) = 1 = u_N+1(t)
> v_0(t) = 3  = v_N+1(t),
> 
> Initial conditio
> u_i(0) = 1+sin (2 pi  x_i  )
> v_i = 3
> For i = 1,?., N
> Time Interval =[t_o, t_end] = [0,10]

Since this is 1+1 dimensional initial value problem, Mathematica can do the 
discretization for you automatically.  Just use

NDSolve[{
     D[u[t, x], t] == 1 - 4 u[t, x] + 0.02 D[u[t, x], x, x] + u[t, x]^3 v[t, x],
     D[v[t, x], t] == 3u[t, x] + 0.02 D[v[t, x], x, x] + u[t, x] ^3v[t, x],
     u[t, 0] == u[t, 1] == 1, u[0, x] == 1 + Sin[2 Pi x],
     v[t, 0] == v[t, 1] == 3, v[0, x] == 3},
   {u, v},
   {t, 0, 10}, {x, 0, 1}]

Mathematica by default uses fourth order differences instead of the second order 
you specified above.  If you really want second order spatial differences with 
the exact spacing you defined, you can use

NDSolve[{D[u[t, x], t] ==
       1 - 4 u[t, x] + 0.02 D[u[t, x], x, x] + u[t, x]^3 v[t, x],
     D[v[t, x], t] == 3u[t, x] + 0.02 D[v[t, x], x, x] + u[t, x] ^3v[t, x],
     u[t, 0] == u[t, 1] == 1, u[0, x] == 1 + Sin[2 Pi x],
     v[t, 0] == v[t, 1] == 3, v[0, x] == 3},
   {u, v},
   {t, 0, 10}, {x, 0, 1}, StartingStepSize -> 1./206, MaxSteps -> {1000, 300},
   DifferenceOrder -> 2]

Interestingly enough, either way you do it, Mathematica is only able to carry 
out the solution out to about t == 0.035 or so.  This appears to be because the 
nonlinearity is causing the solution to form a nonsingularity which appears as a 
cusp with rapidly rising height.

The equations you wrote above would not converge to the solution of a PDE.  This 
is because of the (delta (x))^3  in .02 (u_i-1 ? 2 u_i + u i+1)/(delta (x))^3 
Any finite difference formula for the second derivative on a uniform grid will 
involve (delta (x))^2 -- i.e. squared, not cubed.  At a fixed grid space, the 
extra power will have the effect of making the diffusion coefficient larger by a 
factor of n (which for your range of interest effectively removes to formation 
of the discontinuity).  Presumably the (delta (x))^3 was a typo as was the u_i+1 
instead of v_i+1 in .02 (v_i-1 ? 2 v_i +  u i+1)

However, if you really wanted the cubed power, here is how you could manually do 
the discretization with Mathematica:

n = 205;
X = N[Range[1, n]/(n + 1)];
U[t_] = Map[u[#][t] &, Range[1, n]];
V[t_] = Map[v[#][t] &, Range[1, n]];

eqns = Join[
       Thread[D[U[t], t] ==
           1 - 4 U[t] +
             0.02 ListCorrelate[N[{1, -2, 1} n^3], U[t], {2, 2}, 1] +
             U[t]^3 V[t]],
       Thread[
         D[V[t], t] ==
           3 U[t] + 0.02 ListCorrelate[N[{1, -2, 1} n^3], V[t], {2, 2}, 3] +
             U[t]^3 V[t]],
       Thread[U[0] == 1 + Sin[X]],
       Thread[V[0] == 3 + 0. X]];

NDSolve[eqns, Join[U[t], V[t]], {t, 0, 10}]

> 
> 
> Could i get code in Mathematica (by using Euler of 4 Runge - Kutta..)to 
> solve this ordinary differential equation.

I do not recommend solving this with either an Euler method or a RungeKutta 
method for the reason that the potential formation of discontinuities could make 
the ODEs stiff.  The default method Mathematica uses automatically switches to 
methods appropriate for stiff ODEs when needed.  If you want to use a RungeKutta 
method, you can use

NDSolve[eqns, Join[U[t], V[t]], {t, 0, 10}, Method->RungeKutta]

This uses the Runge-Kutta-Fehlberg 4(5) method.

> I am very happy if you give me help
> With best regards
> Khaled
> 
> 
> _________________________________________________________________
> MSN Photos is the easiest way to share and print your photos: 
> http://photos.msn.com/support/worldwide.aspx
> 
> 



  • Prev by Date: RE: Variable functions
  • Next by Date: How to protect the Mathematica Program
  • Previous by thread: Method of line
  • Next by thread: RE: from Jonathan at MathGroup