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