Блог

Вот относительно простая краевая задача, решаемая с помощью метода съемки и 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)
 

Система сингулярна при r=0 (деление на ноль), и мы позаботились об этом, определив предельную дНк на границе r=0

 dNa = ka/3*Ca
 

Другие методы могут решить эту проблему численно (используя r как небольшое число на границе, деля на r небольшое число).

Решение той же проблемы в Гекко, игнорирующее сингулярную краевую задачу, может быть таким

 #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])
 

Гекко не будет жаловаться на сингулярность на границе и решит эту проблему.

И Питон, и Гекко решат эту задачу, удовлетворяя граничным условиям

Питон

 r=0 => [0.02931475 0.        ]
r=R => [ 0.2        -0.12010739]
 

Гекко

 r=0 => 0.029314856906 0.0
r=R => 0.2 -0.12010738377
 

Я не знаю, как включить сингулярность на границе he в Гекко. С другой стороны, Гекко дал результат, не жаловался на сингулярность и выполнял граничные условия Na(0)=0, Ca(R)=0,2.

Я полагаю, что метод коллокации может успешно избежать проблемы с особенностями на границах, но я хотел бы получить подтверждение того, правильно это или нет в Гекко — просто проигнорировать его.

Что можно было с этим поделать?

С наилучшими пожеланиями, Радован

Один из способов преодолеть большинство проблем с делением на ноль-это изменить уравнение, умножив обе стороны на знаменатель, например:

 #m.Equation(Na.dt() == -2/t*Na-ka*Ca)
m.Equation(t*Na.dt() == -2*Na-t*ka*Ca)
 

Когда t=0 , уравнение таково 0==-2*Na . При начальном условии Na.value=0 с Na = m.Var(Na_0) это уравнение выполняется , даже если Гекко не включает уравнения в начальный момент времени.

Это не проблема, но допустимый диапазон узлов составляет от 2 до 6, поэтому, когда m.options.NODES = 7 фактические используемые узлы равны 6.

Ваша проблема выглядит правильной. Попробуйте m.fix_final(Ca,Ca_R) убедиться, что окончательное условие выполнено.