Here is a relatively simple boundary value problem solved with shooting method and Python
import numpy as np
from scipy.optimize import fsolve
from scipy.integrate import odeint
import matplotlib.pyplot as plt
%matplotlib inline
R = 0.5
ka = 6.4
De = 0.1
Cas = 0.2
def odefun(U, r):
Ca, Na = U
dCa = -Na/De
if np.abs(r) <= 1e-10:
dNa = ka/3*Ca
else:
dNa = -2/r*Na-ka*Ca
return [dCa, dNa]
r = np.linspace(0, R)
Na_0 = 0
Ca_R = 0.2
def objective(x):
U = odeint(odefun, [x, Na_0], r)
u = U[-1,0]-Ca_R
return u
x0 = 0.1 #initial guess
x, = fsolve(objective, x0)
print ("Ca_0=",x)
U = odeint(odefun, [x, Na_0], r)
print ("r=0 =>",U[0])
print ("r=R =>",U[-1])
plt.plot(t.value,Ca.value)
plt.plot(t.value,Na.value)
The system is singular at r=0 (divison by zero) and we took care about this by defining limiting dNa at the boundary r=0
dNa = ka/3*Ca
Other methods are possible to solve this numerically (using r as a small number at the boundary, dividing by r+small number)
Solving the same problem in Gekko ignoring the singular boundary problem might be this way
#Boundary value problem
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
R = 0.5
ka = 6.4
De = 0.1
Cas = 0.2
Na_0 = 0
Ca_R = 0.2
m = GEKKO()
nt = 101
m.time = np.linspace(0,R,nt) # time points
Na = m.Var(Na_0) # known at r=0
Ca = m.Var(fixed_initial=False) # unknown at r=0
pi = np.zeros(nt)
pi[-1]=1
p = m.Param(value=pi)
# create GEEKO equations
t = m.Var(m.time)
m.Equation(t.dt() == 1)
m.Equation(Na.dt() == -2/t*Na-ka*Ca)
m.Equation(Ca.dt() == -Na/De)
m.Minimize(p*(Ca - Ca_R)**2) # minimizing at r=R
# solve ODE
m.options.IMODE = 6
m.options.NODES = 7
m.solve(disp=False)
# plot results
print ("r=0 =>",Ca[0],Na[0])
print ("r=R =>",Ca[-1],Na[-1])
plt.plot(r,U[:,1])
plt.plot(r,U[:,0])
Gekko will not complain about singularity at the boundary and will solve this.
Both Pythone and Gekko will solve this satisfying boundary conditions
Python
r=0 => [0.02931475 0. ]
r=R => [ 0.2 -0.12010739]
Gekko
r=0 => 0.029314856906 0.0
r=R => 0.2 -0.12010738377
I do not know how to include the singularity at he boundary in Gekko. On the other hand, Gekko gave the result, did not complain about the singularity and boundary conditions are satisfied Na(0)=0, Ca(R)=0.2.
I suppose that collocation method can successfully avoid the problem with singularities at the boundaries, but I would like the confirmation whether this is correct or not in Gekko - just to ignore it.
What could be done about this?
Best Regards, Radovan
