MAJOR problem with precision in MAPLE do loops

Ask a Question

Create a Post

Question:MAJOR problem with precision in MAPLE do loops

Posted: Thanos 5 Product: Maple

May 22 2012

Dear Maple Specialists,

I deal with a major MAPLE problem the last months of my PhD.

PROBLEM:
My programm contains a do loop which calculates for given PsiS (purple color) and for a range of Psi1 and Psi2 (blue color) from PsiS to 0, values for the products Qsu1, Qsu2, Q_tot .
The aim to calculate w(RLZ) which is obtained when the blue color loop is completed and Psi2 takes the final value 0.

w(RLZ) depends from Qtot, Qsu1, Qsu2

Qsu1, Qsu2, Q_totare calculated when Psi2 =0 in such a way that always Qtot >= Qsu1, Qsu2.
For PsiS=0.1, Psi1=0.01 and Psi2=0 we get (included calculation)

Qtot = 1.297046409*10^22,
Qsu1 = 1.271370424*10^22
Qsu2 = 1.297046402*10^22

while it is obvious Qtot>Qsu2 MAPLE calculates the difference between then DeltaQ2=Qtot-Qsu2 as NEGATIVE and of the order of 10¹⁴!!!!!!! (DeltaQ2=-4.27e14!!!!!!!!!!)

I conducted the same calculation seperately (see below: seperate calculation) i.e. without do loops and i noticed that the same calculation when is done without loops delivers different values for the IDS,ID1,ID2,IAS,IA1,IA2,ITS,IT1,IT2 and subsequently different values for the sums Qtot, Qsu1, Qsu2 (comparison of the values is found in attached maple sheets)

At the end the differnce DeltaQ2 is correctly calculated by the seperate calculation while my programm -no matter what i tried- calculates Delta Q2 as negative (DeltaQ2 =-4.27e14 !!!!!!!!!) and in any case the Programm calculates different values for all the with red color highlighted elements.

Why maple calculates differently things in the do loop than in the seperate calculation???????
How can i optimise my already written programm to calculate things exactly like the seperate calculation? (My boss proposed mathematica but i dont want to quit MAPLE after one °!"§$%&/ year of good work with it neither learn to use another mathematical software)

I include the programm and the seperate calculation in order to demonstrate precicly my problem.

The elements of the programm and the seperate calculation marked RED are unterlying that there is a disagreement between the two calculation ways (Programm/seperate calculation)
The GREEN elements are calculated correctly from both calculation ways (programm/ seperate calculation)

Thanx in advance for your help

Thanos Tsirimpis

Universität Erlangen Nürnberg @ Angewandte Physik
-----------------------------------------------------------------

PROGRAMM

restart

T := 298:

d := 0.4e-3:

(1)

NC := evalf(0.3e22*T^(3/2));

(2)

Matrix(%id = 171061860)

(3)

(4)

(5)

wiC := Einb/(k*T)-ECB/(k*T);

(6)

Matrix(%id = 150549700)

(7)

for i to 9001 do if `or`(FD32[i, 1] < wiC+uB and wiC+uB < FD32[i+1, 1], FD32[i, 1] = wiC+uB) then iFD32Cb := i else FAIL end if end do;

FD32Cb := FD32[iFD32Cb, 2]+(wiC+uB-FD32[iFD32Cb, 1])*(FD32[iFD32Cb+1, 2]-FD32[iFD32Cb, 2])/(FD32[iFD32Cb+1, 1]-FD32[iFD32Cb, 1]);

(8)

for i to 10001 do if `or`(FD32[i, 1] < wiC+uS and wiC+uS < FD32[i+1, 1], FD32[i, 1] = wiC+uS) then iFD32Cs := i else FAIL end if; if `or`(FD32[i, 1] < wiC+uS and wiC+uS < FD32[i+1, 1], FD32[i, 1] = wiC+uS) then break end if end do:

(9)

for i to 10001 do if `or`(FD32[i, 1] < wiC+u1 and wiC+u1 < FD32[i+1, 1], FD32[i, 1] = wiC+u1) then iFD32C1 := i else FAIL end if; if `or`(FD32[i, 1] < wiC+u1 and wiC+u1 < FD32[i+1, 1], FD32[i, 1] = wiC+u1) then break end if end do:

(10)

IAs := NA*(uB-uS+ln(1/(1+gA*exp(wAi-uS))+gA*exp(wAi-uB)/(1+gA*exp(wAi-uS))));

(11)

ID1 := ND*(u1-uS+ln(1/(1+gD*exp(u1+wiD))+gD*exp(uS+wiD)/(1+gD*exp(u1+wiD)))); IA1 := NA*(u1-uS+ln(1/(1+gA*exp(wAi-uS))+gA*exp(wAi-u1)/(1+gA*exp(wAi-uS)))); IT1 := piecewise(Psi1 > PsiTf, NT*(u1-uS+ln(1/(1+gD*exp(u1+wiT))+gD*exp(uS+wiT)/(1+gD*exp(u1+wiT)))), Psi1 <= PsiTf, NT*(uT-uS+ln(1/(1+gD*exp(uT+wiT))+gD*exp(uS+wiT)/(1+gD*exp(uT+wiT)))))

(12)

ID2 := ND*(u2-uS+ln(1/(1+gD*exp(u2+wiD))+gD*exp(uS+wiD)/(1+gD*exp(u2+wiD)))); IA2 := NA*(u2-uS+ln(1/(1+gA*exp(wAi-uS))+gA*exp(wAi-u2)/(1+gA*exp(wAi-uS)))); IT2 := piecewise(Psi2 > PsiTf, NT*(u2-uS+ln(1/(1+gD*exp(u2+wiT))+gD*exp(uS+wiT)/(1+gD*exp(u2+wiT)))), Psi2 <= PsiTf, NT*(uT-uS+ln(1/(1+gD*exp(uT+wiT))+gD*exp(uS+wiT)/(1+gD*exp(uT+wiT)))))

(13)

Qsu1 := evalf(ID1+IT1-IA1-NC*(FD32C1-FD32Cs));

(14)

invFu1 := sqrt(1/(2*k*T*(Qtot-Qsu1)/(eps4H*epsilon0))):

(15)

SEPERATE CALCULATION

Download Seperate_Calculatio.mw

restart

T := 298:

d := 0.2e-3:

(1)

NC := evalf(0.3e22*T^(3/2)):

wiC := Einb/(k*T)-ECB/(k*T); wVi := EVB/(k*T)-Einb/(k*T); `Φb` := EF/q-Einb/q; uB := q*`Φb`/(k*T); Psi2 := Psi1-0.1e-1; wiD := Einb/(k*T)-ED/(k*T); wAi := EA/(k*T)-Einb/(k*T); wiT := Einb/(k*T)-ET/(k*T); `Φs` := EF/q-Einb/q-PsiS; Phi1 := EF/q-Einb/q-Psi1; Phi2 := EF/q-Einb/q-Psi2; `ΦT` := EF/q-Einb/q-PsiTf; uS := q*`Φs`/(k*T); u1 := q*Phi1/(k*T); u2 := q*Phi2/(k*T); uT := q*`ΦT`/(k*T)

Matrix(%id = 152287568)

(2)

for i to 9001 do if `or`(FD32[i, 1] < wiC+uB and wiC+uB < FD32[i+1, 1], FD32[i, 1] = wiC+uB) then iFD32Cb := i else FAIL end if end do;

FD32Cb := FD32[iFD32Cb, 2]+(wiC+uB-FD32[iFD32Cb, 1])*(FD32[iFD32Cb+1, 2]-FD32[iFD32Cb, 2])/(FD32[iFD32Cb+1, 1]-FD32[iFD32Cb, 1]);

(3)

IDs := evalf(ND*(uB-uS+ln(1/(1+gD*exp(uB+wiD))+gD*exp(uS+wiD)/(1+gD*exp(uB+wiD))))):

ID1 := evalf(ND*(u1-uS+ln(1/(1+gD*exp(u1+wiD))+gD*exp(uS+wiD)/(1+gD*exp(u1+wiD))))); IA1 := evalf(NA*(u1-uS+ln(1/(1+gA*exp(wAi-uS))+gA*exp(wAi-u1)/(1+gA*exp(wAi-uS))))); IT1 := evalf(piecewise(Psi1 > PsiTf, NT*(u1-uS+ln(1/(1+gD*exp(u1+wiT))+gD*exp(uS+wiT)/(1+gD*exp(u1+wiT)))), Psi1 <= PsiTf, NT*(uT-uS+ln(1/(1+gD*exp(uT+wiT))+gD*exp(uS+wiT)/(1+gD*exp(uT+wiT))))))

ID2 := evalf(ND*(u2-uS+ln(1/(1+gD*exp(u2+wiD))+gD*exp(uS+wiD)/(1+gD*exp(u2+wiD))))); IA2 := evalf(NA*(u2-uS+ln(1/(1+gA*exp(wAi-uS))+gA*exp(wAi-u2)/(1+gA*exp(wAi-uS))))); IT2 := evalf(piecewise(Psi2 > PsiTf, NT*(u2-uS+ln(1/(1+gD*exp(u2+wiT))+gD*exp(uS+wiT)/(1+gD*exp(u2+wiT)))), Psi2 <= PsiTf, NT*(uT-uS+ln(1/(1+gD*exp(uT+wiT))+gD*exp(uS+wiT)/(1+gD*exp(uT+wiT))))))