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.
Comments
Post a Comment