奇异常微分方程边值问题

Question

这是一个相对简单的边值问题，用打靶法和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来解决这个问题

dNa = ka/3*Ca

其他方法也可以数值求解（在边界处使用r作为小数值，除以r+小数值）

在Gekko中解决同样的问题忽略奇异边界问题可能是这样的

#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不会抱怨边界的奇异性，而是会解决这个问题

Pythone和Gekko都将解决这个令人满意的边界条件

Python

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

盖柯

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

我不知道如何将Gekko边界的奇点包括在内。另一方面，Gekko给出了结果，没有抱怨奇异性和边界条件满足Na（0）=0，Ca（R）=0.2

我想，配置方法可以成功地避免边界奇点的问题，但我想确认这在Gekko中是否正确——只是忽略它

对此我们能做些什么

致以最良好的祝愿， 拉多万

Answer 1

克服大多数“零除”问题的一种方法是，通过将双方乘以分母来改革等式，例如：

#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)，即使Gekko在初始时间点不包括方程，也满足该方程

这不是问题，但是节点的有效范围是2到6，因此m.options.NODES = 7实际使用的节点是6

你的问题看起来是对的。尝试m.fix_final(Ca,Ca_R)以确保满足最终条件