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