具有奇异性的常微分方程边值问题

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+small数)

在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 = m.Var(Na_0)的初始条件Na.value=0，即使Gekko在初始时间点不包括方程，该方程也是满足的。

这不是问题，但是节点的有效范围是2到6，所以当使用m.options.NODES = 7时，实际使用的节点是6。

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