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