Closed MSJavaScript closed 1 year ago
It is not clear to me that you have demonstrated a shortcoming in SciPy's integrate
module. However, it seems likely that the Python code you provided is not an accurate translation of the Matlab script you also attached. I would encourage you to start again and be more careful as you work your way through the translation. In particular, try to avoid adding in unnecessary functions and classes until you are sure the code you have written works. You may want to consider the below Python snippet as a new starting point:
import numpy as np
from scipy.constants import N_A, e, epsilon_0, h, k
from scipy.integrate import solve_ivp
from scipy.sparse import spdiags
from scipy.sparse.linalg import spsolve
def ion_conc_equ(t, y):
return
T = 298.15
a_M = 3.5e-10
N_M = 1 / np.sqrt(3) / a_M**2
potential = 0.5 * e / k / T
delta_RP = 1.5 * 2.75e-10
epsilon_RP = 6
epsilon_s = 78.5
D = 1e-9
n0 = N_A
Ld = 1e-4
lambda_d = np.sqrt(epsilon_s * epsilon_0 * k * T / e**2 / n0)
tref = lambda_d * Ld / D
Length_limit = Ld / lambda_d
x_cal_1 = np.linspace(0, 1, 70)
x_cal_2 = np.linspace(1.1, 10, 50)
x_cal_3 = np.linspace(11, 100, 50)
x_cal_4 = np.linspace(110, Length_limit, 20)
xmesh = np.hstack((x_cal_1, x_cal_2, x_cal_3, x_cal_4))
nx = xmesh.size
harray = np.diff(xmesh)
diag0 = np.hstack(
(
-1 / harray[0] - lambda_d * epsilon_RP / epsilon_s / delta_RP,
-2 / harray[1:-1] / harray[:-2],
-2 / harray[-1] / harray[-2],
)
)
diagp1 = np.hstack((1 / harray[0], 2 / harray[1:-1] / (harray[1:-1] + harray[:-2])))
diagm1 = np.hstack(
(
2 / harray[0] / (harray[0] + harray[1]),
2 / harray[1:-1] / (harray[1:-1] + harray[2:]),
)
)
A = np.column_stack((np.pad(diagm1, (0, 1)), diag0, np.pad(diagp1, (1, 0))))
# https://stackoverflow.com/a/31900194/10640534
pot_M = spdiags(A.T, [-1, 0, 1], nx - 1, nx - 1).tocsc()
pot_b = np.zeros(nx - 1)
pot_b[0] = -potential * lambda_d * epsilon_RP / delta_RP / epsilon_s
y0 = np.ones((5, nx))
y0[-1, :-1] = spsolve(pot_M, pot_b)
y0[-1, -1] = 0
t_end = 5
tmesh = np.linspace(0.1, t_end, 101)
tmesh = np.hstack(([1e-10, 1e-8, 1e-6, 1e-4, 1e-2], tmesh))
t_cal = tmesh / tref
sol = solve_ivp(
ion_conc_equ, [t_cal[0], t_cal[-1]], y0.flatten(), method="LSODA", t_eval=t_cal
)
It should reproduce the Matlab script until line 74, which is the ode15s
call. For example, you could start by making sure that arrays calculated up until there in both languages are identical. Then, you could translate the system's right-hand side into the ion_conc_equ
function. Once you are sure that the Python code you have written is an accurate translation of the whole Matlab script, then you can assess if the ODE solvers provided by SciPy are behaving well or not.
Thank you for your kindly reply. I complement the python version script basing on your script. It seems that the problem is the jac_sparsity parameter for BDF solver or lband/uband parameter for LSODA solver.
For BDF solver, if I leave jac_sparsity None, the solver can slowly converge to ode15s result. But if I set jac_sparsity, the result never converges (I mean reach the end of time span). For LSODA, the problem is similar, if I set lband/uband, the result also cannot converge. It's weird that ode15s cannot converge too if I set its JPattern parameter. But ode15s still converge much faster that solve_ivp if I don't set the Jacobian matrix.
import numpy as np
#from scipy.constants import N_A, e, epsilon_0, h, k
from scipy.integrate import solve_ivp
from scipy.sparse import diags, spdiags, kron
from scipy.sparse.linalg import spsolve
e = 1.60217663e-19
k = 1.380649e-23
N_A = 6.0221408e+23
h = 6.626e-34
epsilon_0 = 8.8541878128e-12
def current(NHm,NAm,NBm,Phim):
global T, potential, N_M
e_0 = e
k_B = k
E_0 = 0.6
G_00 = 0.4*e_0
eta_HP = (potential - Phim)*k_B*T/e_0 - E_0
jr = e_0*N_M*k_B*T/h*np.exp(-G_00/k_B/T)*(NBm*np.exp(e_0*eta_HP/2/k_B/T) - NAm*NHm*np.exp(-e_0*eta_HP/2/k_B/T))
print(jr)
return jr
def ion_conc_equ(t, y):
global nx, harray, Ld, lambda_d, pot_M
ionsD = np.array([1.0, 1.0, 1.0, 1.0])
ionsZ = np.array([1.0, -1.0, 0.0, 0.0])
ionsNumDen = np.array([1.0, 1.0, 1.0, 1.0])
dydt = np.zeros((5, nx))
data = np.reshape(y, (5, nx), order='F')
ion_c_array = data[0:4, :]
pot_x = data[4, :]
NHm = data[0, 0]
NAm = data[2, 0]
NBm = data[3, 0]
Phim = data[4, 0]
jr = current(NHm,NAm,NBm,Phim)
jr = jr * Ld/e/n0/D
flux = np.array([-jr, 0, -jr, jr])
dphidx = (pot_x[1] - pot_x[0])/harray[0]
c = 0.5*(ion_c_array[:, 0] + ion_c_array[:, 1])
dcdx = (ion_c_array[:, 1] - ion_c_array[:, 0])/harray[0]
FL = ionsD * ionsZ * c * dphidx + ionsD * dcdx
FL = FL* Ld / lambda_d
dydt[0:4, 0] = (FL - flux)/(0.5*harray[0])
for i in range(1, nx-1):
dphidx = (pot_x[i+1] - pot_x[i])/harray[i]
c = 0.5*(ion_c_array[:, i] + ion_c_array[:, i+1])
dcdx = (ion_c_array[:,i+1] - ion_c_array[:, i])/harray[i]
FR = ionsD * ionsZ * c * dphidx + ionsD * dcdx
FR = FR * Ld / lambda_d
dydt[0:4, i] = (FR - FL)/(0.5*(harray[i]+harray[i-1]))
FL = FR
dydt[0:4, nx-1] = ion_c_array[:, nx-1] - ionsNumDen
ion_e_array = dydt[0, :] - dydt[1, :]
ion_e_array = -ion_e_array
pot_b = np.zeros((nx-1, ))
pot_b[0] = 0.25*(ion_e_array[0] + ion_e_array[1])*harray[0]
pot_b[1:] = ion_e_array[1:-1]
dydt[4, 0:-1] = spsolve(pot_M, pot_b)
kk = dydt.flatten(order='F')
#print(kk[122], kk[88], kk[76], kk[31])
#print(kk[9], kk[99], kk[104], kk[199])
#a = input("pause: ")
return kk
T = 298.15
a_M = 3.5e-10
N_M = 1 / np.sqrt(3) / a_M**2
potential = 0.5 * e / k / T
delta_RP = 1.5 * 2.75e-10
epsilon_RP = 6
epsilon_s = 78.5
D = 1e-9
n0 = N_A
Ld = 1e-4
lambda_d = np.sqrt(epsilon_s * epsilon_0 * k * T / e**2 / n0)
tref = lambda_d * Ld / D
Length_limit = Ld / lambda_d
x_cal_1 = np.linspace(0, 1, 70)
x_cal_2 = np.linspace(1.1, 10, 50)
x_cal_3 = np.linspace(11, 100, 50)
x_cal_4 = np.linspace(110, Length_limit, 20)
xmesh = np.hstack((x_cal_1, x_cal_2, x_cal_3, x_cal_4))
nx = xmesh.size
harray = np.diff(xmesh)
diag0 = np.hstack(
(
-1 / harray[0] - lambda_d * epsilon_RP / epsilon_s / delta_RP,
-2 / harray[1:-1] / harray[:-2],
-2 / harray[-1] / harray[-2],
)
)
diagp1 = np.hstack((1 / harray[0], 2 / harray[1:-1] / (harray[1:-1] + harray[:-2])))
diagm1 = np.hstack(
(
2 / harray[0] / (harray[0] + harray[1]),
2 / harray[1:-1] / (harray[1:-1] + harray[2:]),
)
)
A = np.column_stack((np.pad(diagm1, (0, 1)), diag0, np.pad(diagp1, (1, 0))))
# https://stackoverflow.com/a/31900194/10640534
pot_M = spdiags(A.T, [-1, 0, 1], nx - 1, nx - 1).tocsc()
pot_b = np.zeros(nx - 1)
pot_b[0] = -potential * lambda_d * epsilon_RP / delta_RP / epsilon_s
y0 = np.ones((5, nx))
y0[-1, :-1] = spsolve(pot_M, pot_b)
y0[-1, -1] = 0
#y0[0, :] = np.sin(xmesh)
#y0[1, :] = np.cos(0.1*xmesh)
#y0[2, :] = np.sin(np.cos(xmesh+2))
#y0[3, :] = np.exp(np.sin(xmesh))
#y0[4, :] = np.exp(-np.cos(xmesh))
t_end = 5
tmesh = np.linspace(0.1, t_end, 101)
tmesh = np.hstack(([1e-10, 1e-8, 1e-6, 1e-4, 1e-2], tmesh))
t_cal = tmesh / tref
sol = solve_ivp(ion_conc_equ, [t_cal[0], t_cal[-1]], y0.flatten(order='F'), method="LSODA", t_eval=t_cal)
#jac_sparsity = kron(spdiags(np.ones((3, nx)), [-1, 0, 1], nx, nx), np.ones((5, 5)))
#jac_sparsity = None
#sol = solve_ivp(ion_conc_equ, [t_cal[0], t_cal[-1]], y0.flatten(order='F'), method="BDF", t_eval=t_cal, jac_sparsity=jac_sparsity)
Below is the matlab code
clear;
global e_0 k_B T h
e_0 = 1.60217663e-19; % [C] electron charge
k_B = 1.380649e-23; % [J K-1] Boltzmann constant
T = 298.15; % [K] absolute temperature
N_A = 6.0221408e+23; % [mol-1] Avogadro constant
h = 6.626e-34; % [J s] Planck constant
epsilon_0 = 8.8541878128e-12; % [F m-1] vacuum permittivity
% Electrode properties
global N_M potential
a_M = 3.5e-10; % [m] lattice constant of the electrode
N_M = (sqrt(3)*(a_M)^2)^(-1); % [m-2] areal number density of M(111)
potential = 0.5;
% Solution properties
global delta_RP epsilon_RP epsilon_s n0 D BIG BIG2
BIG = 1; % magnify n_A in the calculation to reduce calculation error
BIG2 = 1; % magnify n_B in the calculation to reduce calculation error
delta_RP = 1.5*2.75e-10; % [m] distance from the electrode to the HP
epsilon_RP = 6; % dielectric permittivity of the space between the electrode and the HP
epsilon_s = 78.5; % bulk dielectric permittivity of the water solvent medium
D = 1.00e-9; %[mA2/s] diffusion coefficient
n0 =1*N_A; %[m^-3] number density of total cations (anions) in the solution bulk
% reference values
global Ld lambda
Ld = 100e-6; % [m] distance from the HP to the solution bulk
lambda = sqrt(epsilon_s*epsilon_0*k_B*T/e_0^2/n0); % [m] Debye length
tref= lambda*Ld/D; % [s] time reference
Ld = 1.0e-4;
Length_limit = Ld/lambda;
potential = potential * e_0/k_B/T;
x_cal_1 = linspace(0,1,70); % spatial coordinate of 1st layer
x_cal_2 = linspace(1.1,10,50); % spatial coordinate of 2nd layer
x_cal_3 = linspace(11,100,50); % spatial coordinate of 3rd layer
x_cal_4 = linspace(110,Length_limit,20); % spatial coordinate of 4th layer
xmesh = [x_cal_1, x_cal_2, x_cal_3,x_cal_4]; % assembled spatial coordinate
global nx harray pot_M
nx = length(xmesh);
harray = diff(xmesh);
diag0 = zeros(nx-1, 1);
diagp1 = zeros(nx-2, 1);
diagm1 = zeros(nx-2, 1);
for i = 2:(nx-2)
diag0(i) = -2/(harray(i)*harray(i-1));
diagp1(i) = 2/(harray(i)*(harray(i) + harray(i-1)));
diagm1(i) = 2/(harray(i)*(harray(i) + harray(i+1)));
end
diag0(1) = -(1/harray(1) + lambda*epsilon_RP/epsilon_s/delta_RP);
diag0(nx-1) = -2/(harray(nx-1)*harray(nx-2));
diagp1(1) = 1/harray(1);
diagm1(1) = 2/(harray(1)*(harray(1) + harray(2)));
A = zeros(nx-1,3);
A(1:(end-1),1) = diagm1;
A(:, 2) = diag0;
A(2:end, 3) = diagp1;
pot_M = spdiags(A, [-1, 0, 1], nx-1, nx-1);
pot_b = zeros(nx-1, 1);
pot_b(1) = -potential*lambda*epsilon_RP/(delta_RP*epsilon_s);
y0 = zeros(5, nx);
y0(1:4, :) = 1.0;
y0(5, 1:(end-1)) = pot_M \ pot_b;
% y0(1, :) = sin(xmesh);
% y0(2, :) = cos(0.1*xmesh);
% y0(3, :) = sin(cos(xmesh+2));
% y0(4, :) = exp(sin(xmesh));
% y0(5, :) = exp(-cos(xmesh));
y0 = y0(:);
t_end = 5; % [s]
tmesh = linspace(0.1,t_end,101);
tmesh = [1.0e-10, 1.0e-8, 1.0e-6, 1.0e-4, 1.0e-2, tmesh];
t_cal = tmesh/tref;
S = kron( spdiags(ones(nx,3),-1:1,nx,nx), ones(5,5));
options = odeset('RelTol',1e-6,'AbsTol',1e-9, 'JPattern', S);
options = odeset('RelTol',1e-6,'AbsTol',1e-9);
[t,y] = ode15s(@ion_conc_equ, t_cal, y0, options);
nt = length(t_cal);
NnH = zeros(nt, nx);
Nnn = zeros(nt, nx);
NnA = zeros(nt, nx);
NnB = zeros(nt, nx);
Pphi= zeros(nt, nx);
for i = 1:nt
sol = reshape(y(i, :), 5, nx);
NnH(i, :) = sol(1, :);
Nnn(i, :) = sol(2, :);
NnA(i, :) = sol(3, :);
NnB(i, :) = sol(4, :);
Pphi(i, :) = sol(5, :);
end
f2 = figure(2);
hold on;
for i = 1:100
if(mod(i, 10) == 1)
plot(xmesh*lambda/1e-6, NnA(i,:)*n0/N_A,'LineWidth',2) % 0.1 s
end
end
hold off;
function dydt = ion_conc_equ(t, y)
global nx Ld lambda e_0 n0 D harray pot_M
dydt = zeros(5, nx);
data = reshape(y, 5, nx);
ion_c_array = data(1:4, :);
pot_x = data(5,:);
NHm = data(1, 1);
NAm = data(3, 1);
NBm = data(4, 1);
Phim = data(5, 1);
jr = current(NHm,NAm,NBm,Phim);
jr = jr * Ld/e_0/n0/D;
flux = [-jr, 0, -jr, jr]';
ionsD = [1, 1, 1, 1]';
ionsZ = [1, -1, 0, 0]';
ionsNumDen = [1, 1, 1, 1]';
dphidx = (pot_x(2) - pot_x(1))/harray(1);
c = 0.5*(ion_c_array(:, 1) + ion_c_array(:, 2));
dcdx = (ion_c_array(:, 2) - ion_c_array(:, 1))/harray(1);
FL = ionsD .* ionsZ .* c * dphidx + ionsD.*dcdx;
FL = FL*Ld/lambda;
dydt(1:4, 1) = (FL - flux)/(0.5*harray(1));
for i = 2:(nx-1)
dphidx = (pot_x(i+1) - pot_x(i))/harray(i);
c = 0.5*(ion_c_array(:, i) + ion_c_array(:, i+1));
dcdx = (ion_c_array(:,i+1) - ion_c_array(:, i))/harray(i);
FR = ionsD .* ionsZ .* c .* dphidx + ionsD.*dcdx;
FR = FR*Ld/lambda;
dydt(1:4, i) = (FR - FL)/(0.5*(harray(i)+harray(i-1)));
FL = FR;
end
dydt(1:4, nx) = ion_c_array(:, nx) - ionsNumDen;
ion_e_array = dydt(1, :) - dydt(2, :);
ion_e_array = -ion_e_array;
pot_b = zeros(nx-1, 1);
pot_b(1) = 0.25*(ion_e_array(1) + ion_e_array(2))*harray(1);
pot_b(2:end) = ion_e_array(2:end-1);
dydt(5, 1:(end-1)) = pot_M \ pot_b;
dydt = dydt(:);
% disp([dydt(123), dydt(89), dydt(77), dydt(32)]);
% disp([dydt(10), dydt(100), dydt(105), dydt(200)]);
% pause;
end
function jr = current(NHm,NAm,NBm,Phim)
global T e_0 k_B h N_M potential
E_0 = 0.6; % [V] Equilibrium potential of the reaction at standard state
G_00 = 0.4*e_0; % [J] Activation energy of the reaction at standard equilibrium state
eta_HP = (potential - Phim)*k_B*T/e_0 - E_0; % [V] overpotential
jr = e_0*N_M*k_B*T/h*exp(-G_00/k_B/T).* (NBm.*exp (e_0*eta_HP/2/k_B/T)-NAm.*NHm.*exp(-e_0*eta_HP/2/k_B/T)); % [A2 m-2] current density
disp(jr);
end
To make sure ion_conc_equ function in these two cases produce the same output, I initialisze $y_0$ with some sin or cos functions. Python and matlab outputs are same before 7 or 8 decimal places.
I am not too sure about the Jacobian structure. If you manage to make it work in Matlab, I can help you translate it into Python.
Does the following work better than what you pasted above, or is it still too slow?
from time import perf_counter
import numpy as np
from scipy.constants import N_A, e, epsilon_0, h, k
from scipy.integrate import solve_ivp
from scipy.sparse import spdiags
from scipy.sparse.linalg import spsolve
def current(NHm, NAm, NBm, Phim, E_0=0.6, G_00=0.4 * e):
eta_HP = (potential - Phim) * k * T / e - E_0
return (
(e * N_M * k * T / h)
* np.exp(-G_00 / k / T)
* (
NBm * np.exp(e * eta_HP / 2 / k / T)
- NAm * NHm * np.exp(-e * eta_HP / 2 / k / T)
)
)
def ion_conc_equ(t, y):
print(t)
dydt = np.zeros((5, nx))
data = np.reshape(y, (5, nx), order="F")
ion_c_array = data[:-1]
pot_x = data[-1]
NHm = data[0, 0]
NAm = data[2, 0]
NBm = data[3, 0]
Phim = data[4, 0]
jr = current(NHm, NAm, NBm, Phim)
jr *= Ld / e / n0 / D
flux = np.array([-jr, 0, -jr, jr])
ionsD = np.array([1, 1, 1, 1])
ionsZ = np.array([1, -1, 0, 0])
ionsNumDen = np.array([1, 1, 1, 1])
dphidx = (pot_x[1:] - pot_x[:-1]) / harray
c = (ion_c_array[:, :-1] + ion_c_array[:, 1:]) / 2
dcdx = (ion_c_array[:, 1:] - ion_c_array[:, :-1]) / harray
# dphidx = (pot_x[1] - pot_x[0]) / harray[0]
# c = (ion_c_array[:, 0] + ion_c_array[:, 1]) / 2
# dcdx = (ion_c_array[:, 1] - ion_c_array[:, 0]) / harray[0]
FL = (ionsD * ionsZ * c[:, 0] * dphidx[0] + ionsD * dcdx[:, 0]) * Ld / lambda_d
dydt[:-1, 0] = (FL - flux) * 2 / harray[0]
for i in range(1, nx - 1):
# dphidx = (pot_x[i + 1] - pot_x[i]) / harray[i]
# c = (ion_c_array[:, i] + ion_c_array[:, i + 1]) / 2
# dcdx = (ion_c_array[:, i + 1] - ion_c_array[:, i]) / harray[i]
FR = (ionsD * ionsZ * c[:, i] * dphidx[i] + ionsD * dcdx[:, i]) * Ld / lambda_d
dydt[:-1, i] = (FR - FL) * 2 / (harray[i] + harray[i - 1])
FL = FR
dydt[:-1, -1] = ion_c_array[:, -1] - ionsNumDen
ion_e_array = dydt[1] - dydt[0]
pot_b = np.zeros(nx - 1)
pot_b[0] = (ion_e_array[0] + ion_e_array[1]) * harray[0] / 4
pot_b[1:] = ion_e_array[1:-1]
dydt[-1, :-1] = spsolve(pot_M, pot_b)
return dydt.flatten(order="F")
T = 298.15
a_M = 3.5e-10
N_M = 1 / np.sqrt(3) / a_M**2
potential = 0.5 * e / k / T
delta_RP = 1.5 * 2.75e-10
epsilon_RP = 6
epsilon_s = 78.5
D = 1e-9
n0 = N_A
Ld = 1e-4
lambda_d = np.sqrt(epsilon_s * epsilon_0 * k * T / e**2 / n0)
tref = lambda_d * Ld / D
Length_limit = Ld / lambda_d
x_cal_1 = np.linspace(0, 1, 70)
x_cal_2 = np.linspace(1.1, 10, 50)
x_cal_3 = np.linspace(11, 100, 50)
x_cal_4 = np.linspace(110, Length_limit, 20)
xmesh = np.hstack((x_cal_1, x_cal_2, x_cal_3, x_cal_4))
nx = xmesh.size
harray = np.diff(xmesh)
diag0 = np.hstack(
(
-1 / harray[0] - lambda_d * epsilon_RP / epsilon_s / delta_RP,
-2 / harray[1:-1] / harray[:-2],
-2 / harray[-1] / harray[-2],
)
)
diagp1 = np.hstack((1 / harray[0], 2 / harray[1:-1] / (harray[1:-1] + harray[:-2])))
diagm1 = np.hstack(
(
2 / harray[0] / (harray[0] + harray[1]),
2 / harray[1:-1] / (harray[1:-1] + harray[2:]),
)
)
A = np.column_stack((np.pad(diagm1, (0, 1)), diag0, np.pad(diagp1, (1, 0))))
# https://stackoverflow.com/a/31900194/10640534
pot_M = spdiags(A.T, [-1, 0, 1], nx - 1, nx - 1).tocsc()
pot_b = np.zeros(nx - 1)
pot_b[0] = -potential * lambda_d * epsilon_RP / delta_RP / epsilon_s
y0 = np.ones((5, nx))
y0[-1, :-1] = spsolve(pot_M, pot_b)
y0[-1, -1] = 0
t_end = 5
tmesh = np.linspace(0.1, t_end, 101)
tmesh = np.hstack(([1e-10, 1e-8, 1e-6, 1e-4, 1e-2], tmesh))
t_cal = tmesh / tref
begin = perf_counter()
sol = solve_ivp(
ion_conc_equ,
[t_cal[0], t_cal[-1]],
y0.flatten(order="F"),
method="LSODA",
t_eval=t_cal,
)
print(perf_counter() - begin)
print(sol.t, sol.y)
@Patol75 Thanks for your help. I think I know where I am wrong. In fact I write my code basing on a matlab function pdepe. This function converts parabolic and elliptic PDEs to ODEs and solve them by ode15s, If an elliptic PDE exists (for example Poisson equation here), the eqution becomes algebraic differential equation. ODE15s can decouple the equations somehow (I don't know the details), so pdepe function can set a Jacobian pattern like But I my function, because I solve the poisson equation manually with spsolve, the Jacobian is not sparse matrix anymore. So if I set Jacobian matrix with the form above, it will not work.
Describe your issue.
I am trying to solve pne dimensional Possion-Nernst-Planck (PNP) Equation with solve_ivp, it looks like Here, $\psi$ is electrode potential, $c_i$ are concentrations of ions. Here is an example in matlab. matlab_cc.txt The function current calculate the currenr density of a reaction. The result will be used to calculate the Robin condition of left boundary. It's a little complex, but jr should decrease during the integration of equation. And ode15s can reach the end of the time span I want to calculate.
I also write a Python version and use solve_ivp to solve this equation. mkm_electro4.py.txt However, the current density jr decrease slowly. If you print the t variable in ions_conc_equ, you will find the time step is very small, t can only reach 8.0e-7. I tried BDF, LSODA solver, but both of them can not reach the end of the time span I want to calculate.
In short, PNP equation is a stiff equation, but BDF or LSODA solver fail to solve it.
Reproducing Code Example
Error message
SciPy/NumPy/Python version information
1.7.1 1.21.2 sys.version_info(major=3, minor=8, micro=10, releaselevel='final', serial=0)