Solving systems of differential equations
in [Mathematica]
|
From: Yun Zhao on 24 Nov 2009 05:48 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?
From: Bob Hanlon on 25 Nov 2009 02:34 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(a)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?
From: dh on 25 Nov 2009 02:35 Yun Zhao 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? > > Hi, I can not see a difference: rho0 = 0.0258; d0 = 0.0123; p0 = 30000; n0 = 0; system = {p'[t] == -(rho + d)*p[t], n'[t] == rho*p[t] + (rho - d)*n[t]}; sol = DSolve[system, {p, n}, t][[1]] /. {rho -> rho0, d -> d0, C[1] -> p0, C[2] -> n0} Plot[Evaluate[{p[t], n[t]} /. sol], {t, 0, 120}] sol4 = NDSolve[{p'[t] == -(rho + d)*p[t], n'[t] == rho*p[t] + (rho - d)*n[t], p[0] == p0, n[0] == n0} /. {rho -> rho0, d -> d0}, {p, n}, {t, 0, 120}] Plot[Evaluate[{p[t], n[t]} /. sol4], {t, 0, 120}] Daniel
|
Pages: 1 Prev: How to get the output of FindFit from MathLink in Next: Bug in Plot3D? |