MathGroup Archive 2009

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

Search the Archive

Re: Solving systems of differential equations

  • To: mathgroup at smc.vnet.net
  • Subject: [mg105242] Re: [mg105197] Solving systems of differential equations
  • From: Bob Hanlon <hanlonr at cox.net>
  • Date: Wed, 25 Nov 2009 02:33:57 -0500 (EST)
  • Reply-to: hanlonr at cox.net

You didn't show all of your code. The following gives good results for me

system = {
   p'[t] == -(rho + d)*p[t],
   n'[t] == rho*p[t] + (rho - d)*n[t],
   n[0] == 0, p[0] == p0};

sub = {rho -> 258/10000, d -> 123/10000, p0 -> 30000};

sol = DSolve[system, {p, n}, t][[1]];

Plot[{p[t], n[t]} /. sol /. sub, {t, 0, 120}]

sol4 = NDSolve[system /. sub,
    {p[t], n[t]}, {t, 0, 120}][[1]];

Plot[{p[t], n[t]} /. sol4, {t, 0, 120}]


Bob Hanlon

---- Yun Zhao <yun.m.zhao at gmail.com> wrote: 

=============
I have a question as to an apparent inconsistency I noticed while using
Mathematica to solve differential equations, I have 2 ODE like this, and I
solve them symbolically using Mathematica 7.0:

system = {p'[t] == -(rho + d)*p[t],
n'[t] == rho*p[t] + (rho - d)*n[t]};
sol = DSolve[system, {p, n}, t][[1]]

I get the following solution:

{p -> Function[{t}, E^((-d - rho) t) C[1]],
n -> Function[{t},
1/2 (-E^((-d - rho) t) + E^((-d + rho) t)) C[1] +
E^((-d + rho) t) C[2]]}

It looks to be correct.

Then, manually, I know p(0)=C1, n(0)=0, so C1=p0, C2=0.

I know the values for rho and d. So for t from 0 to 120, I can calculate and
plot the values for p(t) and n(t).

Up to now, this is all good.

Now I want to solve the system of ODEs numerically (even though I have
already solved it, I have good reasons to want to learn to solve it
numerically, as I am going to add more complexity to this model), so I did
the following:

sol4 = NDSolve[{p'[t] == -(0.0258 + 0.0123)*p[t],
n'[t] == 0.0258*p[t] + (0.0258 - 0.0123)*n[t], p[0] == 30000,
n[0] == 0}, {p, n}, {t, 0, 120}]

Now I get the following solution:

{{p-> InterpolatingFunction[{{0.`,120.`}},"<>"],
n->InterpolatingFunction[{{0.`,120.`}},"<>"]}}

It may be correct, but I then plotted it using

Plot[Evaluate[{n[t]} /. First[sol4]], {t, 0, 120},
PlotRange -> All]

Now my question is:

I generated plot #1 of n(t) by solving the system of ODEs first
symbolically, then manually inputting in values (C1, C2, rho, d).
I generated plot #2 of n(t) by solving the system of ODEs numerically, then
had mathematica plot the solutions of and n(t) for me.

Why does plot #1 and #2 look so different? By that I mean, n(t) in plot 2 is
much lower than n(t) in plot 1.

Could you please explain what I did wrong? Thank you.

I know mathgroup doesn't like attachment, but is there anyway I can show
these plots?



  • Prev by Date: Re: Delay Differential Equations
  • Next by Date: Re: Solving systems of differential equations
  • Previous by thread: Solving systems of differential equations
  • Next by thread: Re: Solving systems of differential equations