I am badly stuck with this program I am writing. I want output of program below to be precise up to 25 digits after decimal point. I am supplying all inputs with precision 25.
But after each iteration of main for loop, precision keeps on decreasing by some points. And as a result in about 100 iterations precision comes down to 0.
I tried everything that I can think of including explicitly setting precision using SetPrecision
in f1,f2,f3,N1,N2,N3,a
. That would produce output but it's not correct as SetPrecision
is just padding the digits from right.
Appreciate any help/suggestions. Thankx.
Input:
Clear["Global`*"]; order = 0.89`25; parameter = 27.3`25; ic = {1000000001/100000000`25.,
10, 10}; SIZE = 100;
Program:
AbsoluteTiming[alph = bet = gam = order; h = 2/100; omeg = -2667/1000;
mu = 10; A = parameter; B = 1;
Array[y1, SIZE, 0]; y1[0] = ic[[1]];
Array[y2, SIZE, 0]; y2[0] = ic[[2]]; Array[y3, SIZE, 0];
y3[0] = ic[[3]];
a[0 _, k_] := (k - 1)^(alph + 1) - (k - 1 - alph) k^alph;
a[j_, k_] := (k - j + 1)^(alph + 1) + (k - 1 - j)^(alph + 1) -
2 (k - j)^(alph + 1) /; (1 <= j <= k - 1);
a[j_, k_] := 1 /; j == k;
l1 := h^alph/Gamma[alph + 2]; l2 := h^bet/Gamma[bet + 2];
l3 := h^gam/Gamma[gam + 2];
f1[t_, n_] :=
l1*Sum[a[j, n]*omeg*y1[j], {j, 0, n - 1}] -
l1*Sum[a[j, n]*y2[j]*y2[j], {j, 0, n - 1}];
f2[t_, n_] :=
mu + l2*Sum[a[j, n]*mu*y3[j], {j, 0, n - 1}] -
l2*Sum[a[j, n]*mu*y2[j], {j, 0, n - 1}];
f3[t_, n_] :=
mu + l3*Sum[a[j, n]*A*y2[j], {j, 0, n - 1}] -
l3*Sum[a[j, n]*B*y3[j], {j, 0, n - 1}] +
l3*Sum[a[j, n]*y1[j]*y2[j], {j, 0, n - 1}];
N1[u1_, u2_, u3_] := l1*(omeg*u1 - u2*u2);
N2[u1_, u2_, u3_] := l2*(mu*u3 - mu*u2);
N3[u1_, u2_, u3_] := l3*(A*u2 - B*u3 + u1*u2);
For[i = 1, i <= SIZE - 1, i++, (*main loop*)
y10 = f1[h*i, i];
y20 = f2[h*i, i];
y30 = f3[h*i, i];
y11 = N1[y10, y20, y30];
y21 = N2[y10, y20, y30];
y31 = N3[y10, y20, y30];
y12 = N1[y10 + y11, y20 + y21, y30 + y31] - N1[y10, y20, y30];
y22 = N2[y10 + y11, y20 + y21, y30 + y31] - N2[y10, y20, y30];
y32 = N3[y10 + y11, y20 + y21, y30 + y31] - N3[y10, y20, y30];
y1[i] = y10 + y11 + y12;
y2[i] = y20 + y21 + y22;
y3[i] = y30 + y31 + y32;];
xx = Table[y1[i], {i, 0, SIZE - 1}];
yy = Table[y2[i], {i, 0, SIZE - 1}];
zz = Table[y3[i], {i, 0, SIZE - 1}];
]
Output error: General::ovfl: Overflow occurred in computation. >>
Further it tells me there are no significant digits remaining to display.
Output:
In[78]:= xx
Out[78]= {1.000000000000000000000000, -3.35160956937173040383529, \
-7.6063107290506868684570, -13.965998348176594403068, \
-23.217734641185881084877, -34.75517400137475702813, \
-45.6273435387493774249, -52.0072035922437294633, \
-52.794258502857519941, -49.988719962104854383, \
-45.92876611511978598, -41.92397139618361145, -38.44386042357929135, \
-35.5347038434195767, -33.090746410617466, -30.988018715710107, \
-29.13155464170869, -27.46064718989206, -25.9403613009447, \
-24.5520756248237, -23.286816115939, -22.14149256760, \
-21.11722129330, -20.2189730943, -19.4560502901, -18.843109037, \
-18.401537933, -18.16095921, -18.16039845, -18.4482145, -19.0791605, \
-20.106105, -21.563607, -23.44219, -25.65784, -28.0311, -30.298, \
-32.166, -33.40, -33.90, -33.7, -32.9, -32., -30., -3.*10^1, \
-3.*10^1, -0.*10^1, 0.*10^2, 0.*10^3, 0.*10^7, 0.*10^19, 0.*10^63,
0.*10^219, 0.*10^766, 0.*10^2695, 0.*10^9485, 0.*10^33394,
0.*10^117589, 0.*10^414071, 0.*10^1458094, 0.*10^5134494,
0.*10^18080495, 0.*10^63668287, 0.*10^224200219, Overflow[],
Overflow[], Overflow[], Overflow[], Overflow[], Overflow[],
Overflow[], Overflow[], Overflow[], Overflow[], Overflow[],
Overflow[], Overflow[], Overflow[], Overflow[], Overflow[],
Overflow[], Overflow[], Overflow[], Overflow[], Overflow[],
Overflow[], Overflow[], Overflow[], Overflow[], Overflow[],
Overflow[], Overflow[], Overflow[], Overflow[], Overflow[],
Overflow[], Overflow[], Overflow[], Overflow[], Overflow[]}
Precision:
In[79]:= Precision[xx]
Out[79]= 0.
Program works fine with MachinePrecision
but that is not output I want.
$MinPrecission
to desired value? Like$MinPrecission=25
?? $\endgroup$