具有奇点的 ODE 边值问题

问题描述

这里是一个比较简单的边值问题，用射击方法和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,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,0])

Gekko 不会抱怨边界处的奇点，并会解决这个问题。

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

蟒蛇

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 这是否正确 - 只是忽略它。

有什么办法可以解决这个问题？

最好的问候， 拉多万

解决方法

克服大多数被零除问题的一种方法是通过将两边都乘以分母来改进等式，例如：

#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) 以确保满足最终条件。