具有奇点的 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 这是否正确 - 只是忽略它。

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

最好的问候，拉多万

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

gekko python singular