Bonjour et bonne année à tous !

J'ai récemment tenté de convertir en python un programme créé par un collègue en C++ (ci-après le début) afin de le modifier un peu et de le rendre plus automatique.
Code C++ : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
void pause();
 
 
int main(int argc, char *argv[])
 
{
   int Nmax = 4000;
   int I = 0;
   int J = 0;
   int N = 0;
   long double Psif[8100];
   long double PsiTotal[8100];
   long double Rqd[8100];
   long double Pota[8100];
   long double Dtau;
   long double Csi;
   long double Tau = 0;
   long double Tauf = 0;
   long double u;
   long double Pi = 3.14159265358979323846264338327950288419716939937510;
   long double Sum;
   long double Ct1;
   long double Ct2;
   long double Ct3;
   long double Ct4;
   long double Ct5;
   long double Lambda;
   long double Alpha = 0.5;
   long double PsiProv = 0;
   long double CsiPrime;
   long double Beta = 0; // Zero pour lest tests
   long double Gamma = 0; // Zero pour les tests
   long double Rau;
   long double CsiP;
   long double PsiP;
   long double CsiM;
   long double PsiM;
   double Ei=-1.0;
   double Ef=1.0;
   double E0=0.39;
   double S=0.0706858;
   S=S/10000;
   double C0=0.001;
   C0=C0*1000;
   double T=294.65;
   double D=0.0001;
   D=D/10000;
   double Ru=165;
   double Cd=0.000017;
   double ks=20;
   ks=ks/100;
   double v=0.05;
   long double FRT;
   long double FSC0;
   long double DFRT;
   long double F = 96485.3329;
   long double R = 8.3144621;
   long double Ep1;
   long double Ep2;
   long double Ip1;
   long double Ip2;
    FRT = F / R / T;
    FSC0 = F * S * C0;
    DFRT = D * FRT;
 
    //parameters are defined for oxidation
    u = -(Ei - E0) * FRT;
    CsiM = (Ef - E0) * FRT;
    Lambda = ks / sqrt (v * DFRT);
    Rau = FRT * Ru * FSC0 * sqrt (v *DFRT);
    Gamma = Cd * sqrt (v) / FSC0 /sqrt (DFRT);
 
    printf (" Ei = %f ; Ef = %f ; S = %e ;  C0 = %e  \n", Ei, Ef, S, C0);
    printf (" T = %f ; D = %e ; Ru = %e ;  Cd = %e  \n", T, D, Ru, Cd);
    printf (" ks = %f ; v = %f \n", ks, v);
 
 
    Psif[0] = 0;
    Psif[1] = 0;
    Dtau = (CsiM + u) / Nmax; // calculation of the time step.
    Tau = Dtau;
    Csi = -u;
    Beta = Rau * Gamma;
    Ct3 = Beta/Dtau;
 
    printf (" Lambda = %LF ; Rau = %Lf ; Gamma = %Lf ; Dtau = %Lf \n", Lambda, Rau, Gamma, Dtau);
 
    for (I = 1 ; I <= 2.0 * Nmax + 1; I++)
 
        {
            Rqd[I] = (powl ((I + 1.0) , 1.5)) - (2.0 * powl (I, 1.5)) + (powl ((I - 1) , 1.5));
        }
 
 // Forward Scan
 
 
    for (I = 1 ; I <= Nmax ; I++)
        {
            Ct1 = expl (-Alpha * Csi) / Lambda;
            Ct2 = 4.0 / 3.0 * sqrtl (Dtau / Pi) * (1.0 + expl (-Csi));
            Pota[I] = Csi;
 
            Sum = 0;
 
            for (N = 1 ; N < I ; N++)
 
                {
                    Sum = Sum + (Psif[N] * (1.0 + Ct3) - Psif[N - 1] * Ct3) * Rqd[I - N];
                }
 
            N = 0;
 
            Psif[I] = Psif[I - 1];
 
            do
 
                {
                    Psif[I] = (Psif[I] + 100.0 * PsiProv) /101.0; // Il faut calmer la convergence
                    PsiProv = Psif[I];
                    CsiPrime = Csi - Rau * Psif[I] - Beta * (1.0 - expl (-Tau/Beta));
                    Ct1 = expl (-Alpha * CsiPrime) / Lambda;
                    Ct2 = 4.0 / 3.0 * sqrtl (Dtau / Pi) * (1.0 + expl (-CsiPrime));
 
                    Psif[I] = (1.0 - Sum * Ct2 + (Ct1 * Ct3 + Ct2 * Ct3) * Psif[I - 1] ) / (Ct1 + Ct2 + Ct1 * Ct3 + Ct2 * Ct3);
 
                    N++;
 
                }
 
            while ( fabs ( Psif[I] - PsiProv ) > 0.0001 );
 
            PsiTotal[I] = Psif[I] + Gamma * (1.0 - expl (-Tau / Beta));
 
            if ((I % 200) == 0)
 
            {
 
                printf (" I= %d Csi= %Lf Psif= %Lf PsiTotal= %Lf N= %d \n", I, Csi, Psif[I], PsiTotal[I], N);
 
            }
 
            Csi = Csi + Dtau;
            Tau = Tau + Dtau;
        }
 
    Tauf = Tau - Dtau;
}

Cependant, je me heurte à un problème de boucle infinie dans mon programme python et qui n'existe pas en C++.
Code Python : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
import math
import numpy as np
import statistics
 
Nmax = int(8000)
I = int(0)
J = int(0)
N = int(0)
 
Psif = [0,0]
PsiTotal = [0,0]
Rqd = [0,0]
Pota = [0,0]
for i in range(3*Nmax):
    Psif.append(np.longdouble(0e0))
    PsiTotal.append(np.longdouble(0e0))
    Rqd.append(np.longdouble(0e0))
    Pota.append(np.longdouble(0e0))
 
Tau = float(0)
Tauf = 0
PsiProv = np.longdouble(0e0)
F = float(96485.3329)
R = float(8.3144621)
 
Ei=-1
Ef=1
E0=0.39
S=0.0706858
S=S/10000
C0=0.001
C0=C0*1000
T=294.65
D=0.0001
D=D/10000
Ru=165
Cd=0.000017
ks=20
ks=ks/100
Alpha = 0.5
v=0.05
 
FRT=F/R/T
FSC0=F*S*C0
DFRT=D*FRT
u=-(Ei-E0)*FRT
CsiM=(Ef-E0)*FRT
Lambda=ks/math.sqrt(v*DFRT)
Rau=FRT*Ru*FSC0*math.sqrt(v*DFRT)
Gamma=Cd*math.sqrt(v)/FSC0/math.sqrt(DFRT)
 
"calculation of step time"
Dtau = (CsiM+u)/Nmax
Tau=Dtau
Csi=-u
Beta=Rau*Gamma
Ct3=Beta/Dtau
 
print(f"Lambda = {Lambda} ; Rau = {Rau} ; Gamma = {Gamma} ; Dtau = {Dtau}")
 
for i in range(2*Nmax+1):
    Rqd[i]=math.pow((i+2),1.5)-(2*math.pow(i+1,1.5))+(math.pow((i),1.5))
 
"Forward Scan"
 
for i in range(Nmax+1):
    if i==0:
        pass
    else:
        Ct1=math.exp(-Alpha*Csi/Lambda)
        Ct2=np.longdouble(4/3*math.sqrt(Dtau/math.pi)*(1+math.exp(-Csi)))
        Pota[i]=Csi
        Sum=0
        for N in range(i):
            if N==0:
                pass
            else :
                Sum=Sum+(Psif[N]*(1+Ct3)-Psif[N-1]*Ct3)*Rqd[i-N]
        N=int(0)
        Psif[i]=Psif[i-1]
        while True:
            Psif
            Psif[i] = np.longdouble((Psif[i] + 100*PsiProv)/101)
            PsifProv = np.longdouble(Psif[i])
            CsiPrime = Csi - Rau*Psif[i] -Beta*(1+math.exp(-Tau/Beta))
            Ct1 = math.exp(-Alpha*CsiPrime)/Lambda
            Ct2 = np.longdouble(4/3*math.sqrt(Dtau/math.pi)*(1+math.exp(-CsiPrime)))
            Psif[i] = np.longdouble((1-Sum*Ct2+(Ct1*Ct3+Ct2*Ct3)*Psif[i-1])/(Ct1+Ct2+Ct1*Ct3+Ct2*Ct3))
            N=N+1
            if abs(Psif[i]-PsiProv)<=0.0001:
                break
        PsiTotal[i]=Psif[i]+Gamma*(1-math.exp(-Tau/Beta))
        if i%200==0 and i!=0:
            print(f"I= {i}, Csi= {Csi}, Psif= {Psif[i]}, PsiTotal= {PsiTotal[i]}, N= {N}")
        Csi=Csi+Dtau
        Tau=Tau+Dtau
Tauf=Tau-Dtau

Quelqu'un aurait la solution svp ?