MathGroup Archive 2011

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

Search the Archive

problem of numerical error while using NDsolve

  • To: mathgroup at smc.vnet.net
  • Subject: [mg118316] problem of numerical error while using NDsolve
  • From: tarun dutta <tarunduttaz at gmail.com>
  • Date: Thu, 21 Apr 2011 03:12:17 -0400 (EDT)

see the following program..

n = 5;
p = 89/1000;
q = 5/10;
z = 2;
a[-1][t_] = 0; b[-1][t_] = 0; g[-1][t_] = 0; h[-1][t_] = 0;
a[n + 1][t_] = 0; b[n + 1][t_] = 0; g[n + 1][t_] = 0; h[n + 1][t_] = 0;
eqn1 = Table[-b[i]'[
      t] == ((1/2)*i (i - 1) - q*i) a[i][
       t] - (z*p)*(Sqrt[i]*
         a[i - 1][
          t] (Sum[Sqrt[
             i]*(g[i - 1][t]*g[i][t] + h[i - 1][t]*h[i][t]), {i, 0,
            n}]) - Sqrt[i]*
         b[i - 1][
          t] (Sum[Sqrt[
             i]*(g[i - 1][t]*h[i][t] - h[i - 1][t]*g[i][t]), {i, 0,
            n}]) + Sqrt[i + 1]*
         a[i + 1][
          t]*(Sum[Sqrt[
             i + 1]*(g[i + 1][t]*g[i][t] + h[i + 1][t]*h[i][t]), {i,
            0, n}]) -
        Sqrt[i + 1]*
         b[i + 1][
          t] (Sum[Sqrt[
             i + 1] (g[i + 1][t]*h[i][t] - g[i][t]*h[i + 1][t]), {i,
            0, n}])), {i, 0, n}];
eqn2 = Table[
   a[i]'[t] == ((1/2)*i (i - 1) - q*i) b[i][
       t] - (z*p)*(Sqrt[i]*
         b[i - 1][
          t]*(Sum[Sqrt[
             i]*(g[i - 1][t]*g[i][t] + h[i - 1][t]*h[i][t]), {i, 0,
            n}]) + Sqrt[i]*
         a[i - 1][
          t] (Sum[Sqrt[
             i]*(g[i - 1][t]*h[i][t] - h[i - 1][t]*g[i][t]), {i, 0,
            n}]) + Sqrt[i + 1]*
         b[i + 1][

          t]*(Sum[Sqrt[
             i + 1]*(g[i + 1][t]*g[i][t] + h[i + 1][t]*h[i][t]), {i,
            0, n}]) +
        Sqrt[i + 1]*
         a[i + 1][
          t] (Sum[Sqrt[
             i + 1] (g[i + 1][t]*h[i][t] - g[i][t]*h[i + 1][t]), {i,
            0, n}])), {i, 0, n}];
eqn3 = Table[-h[i]'[
      t] == ((1/2)*i (i - 1) - q*i) g[i][
       t] - (z*p)*(Sqrt[i]*
         g[i - 1][
          t] (Sum[Sqrt[
             i]*(a[i - 1][t]*a[i][t] + b[i - 1][t]*b[i][t]), {i, 0,
            n}]) - Sqrt[i]*
         h[i - 1][
          t] (Sum[Sqrt[
             i]*(a[i - 1][t]*b[i][t] - b[i - 1][t]*a[i][t]), {i, 0,
            n}]) + Sqrt[i + 1]*
         g[i + 1][
          t]*(Sum[Sqrt[
             i + 1]*(a[i + 1][t]*a[i][t] + b[i + 1][t]*b[i][t]), {i,
            0, n}]) -
        Sqrt[i + 1]*
         h[i + 1][
          t] (Sum[Sqrt[
             i + 1] (a[i + 1][t]*b[i][t] - a[i][t]*b[i + 1][t]), {i,
            0, n}])), {i, 0, n}];
eqn4 = Table[
   g[i]'[t] == ((1/2)*i (i - 1) - q*i) h[i][
       t] - (z*p)*(Sqrt[i]*
         h[i - 1][
          t]*(Sum[Sqrt[
             i]*(a[i - 1][t]*a[i][t] + b[i - 1][t]*b[i][t]), {i, 0,
            n}]) + Sqrt[i]*
         g[i - 1][
          t] (Sum[Sqrt[
             i]*(a[i - 1][t]*b[i][t] - b[i - 1][t]*a[i][t]), {i, 0,
            n}]) +
        Sqrt[i + 1]*
         h[i + 1][
          t]*(Sum[Sqrt[
             i + 1]*(a[i + 1][t]*a[i][t] + b[i + 1][t]*b[i][t]), {i,
            0, n}]) +
        Sqrt[i + 1]*
         g[i + 1][
          t] (Sum[Sqrt[
             i + 1] (a[i + 1][t]*b[i][t] - a[i][t]*b[i + 1][t]), {i,
            0, n}])), {i, 0, n}];
eqns = Flatten[Join[eqn1, eqn2, eqn3, eqn4]];
bcs = {a[0][0] == -0.039598164614752206,
   a[1][0] == 0.9140241822877686, a[2][0] == -0.04245064788789195,
   a[3][0] == 0.0005070790789081469,
   a[4][0] == -2.01234547273147*10^-6,
   a[5][0] == 4.798752010776119*10^-8,
   b[0][0] == -0.004220040150896003, b[1][0] == 0.39974207839230497,
   b[2][0] == -0.03711109059923467, b[3][0] == 0.0008340810545908181,
   b[4][0] == -8.46285047807362*10^-6,
   b[5][0] == 8.31744477716169*10^-8, g[0][0] == 0.03139062493673353,
   g[1][0] == 0.5648357254746523, g[2][0] == 0.016433486856070365,
   g[3][0] == -0.000010248235328447107,
   g[4][0] == -2.686298302654745*10^-6,
   g[5][0] == -9.140037688766517*10^-8,
   h[0][0] == 0.02450346589620122, h[1][0] == 0.8223105463822804,
   h[2][0] == 0.05393670244886415, h[3][0] == 0.0009755470396067131,
   h[4][0] == 8.195603941693097*10^-6,
   h[5][0] == 3.8526099618300414*10^-8};
var = Join[Table[a[i], {i, 0, n}], Table[b[i], {i, 0, n}],
   Table[g[i], {i, 0, n}], Table[h[i], {i, 0, n}]];
sol = NDSolve[Flatten[{eqns, bcs}], var, {t, 0, 1500},
   MaxSteps -> 10000000, AccuracyGoal -> 8, PrecisionGoal -> 8];
delo = z*Abs[
    Sum[Sqrt[
       i]*(g[i][t]*g[i - 1][t] + h[i - 1][t]*h[i][t] +
        I*(h[i][t]*g[i - 1][t] - h[i - 1][t]*g[i][t])), {i, 0, n}]];
Plot[Evaluate[delo /. sol], {t, 0, 1500}, PlotRange -> All]
argo = z*Arg[
    Sum[Sqrt[
       i]*(g[i][t]*g[i - 1][t] + h[i - 1][t]*h[i][t] +
        I*(h[i][t]*g[i - 1][t] - h[i - 1][t]*g[i][t])), {i, 0, n}]];
arge = z*Arg[
    Sum[Sqrt[
       i]*((a[i][t]*a[i - 1][t] + b[i - 1][t]*b[i][t]) +
        I*(b[i][t]*a[i - 1][t] - b[i - 1][t]*a[i][t])), {i, 0, n}]];
argf = Abs[{argo - arge}];
Plot[Evaluate[argf /. sol], {t, 0, 1500}, PlotRange -> All]

now have a look on the plots(the output of this program).there are so
much disturbance( figure in the plot is not smooth like sine
wave)which is evolved after T=100,.similar disturbance is for second
plot after t=100,there are so much peaks.I think this is due to
numerical error.I have already used various 'starting stepsize' and
changing the 'precisiongoal' and 'accuracygoal' in ''ndsolve' line but
failed to get expected thing.
is there any way of removing this kind of numerical error?kindly give
a valuable insight.
thanks in advance
tarun


  • Prev by Date: Re: trouble printing to PDF
  • Next by Date: Re: Converting list of subscripted variables into an array
  • Previous by thread: RectangleWave[ ] Notebook containing the definitions for who ever will need it
  • Next by thread: Mathematica code for fractal dimension