Python有条件地求解时滞微分方程

Question

我使用dde23 of pydelay软件包来求解一个时滞微分方程。我的问题是:如何有条件地编写方程？例如，目标方程有两个选项：

when x>1, dx/dt=0.25 * x(t-tau) / (1.0 + pow(x(t-tau),10.0)) -0.1*x
otherwise, dx/dt=0.25 * x

我尝试了两种方法，但似乎两种方法都没有用：

1. 方法1没有抱怨，但它的if语句没有被解释。
2. 方法2产生了以下错误： 在函数'double f中找到可执行的c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp：( c:\mingw\bin\g++.exe )c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:734:错误：'x‘不能用作函数c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:734:错误：'x’不能用作函数'PyObject* compiled_func(PyObject*，c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:878:警告:有符号整数表达式和无符号整数表达式的比较c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:883:警告:有符号整数表达式和无符号整数表达式的比较c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:774:警告:未使用的变量'Nhistx_ar‘c：\docume~1\thao\局部变量~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:775:警告:未使用变量'Shistx_ar‘c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:776:警告:未使用变量'Dhistx_ar’c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:785:警告:未使用变量'NVhistx_ar‘c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:786:警告:未使用变量'SVhistx’警告c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:787:警告:未使用变量'DVhistx_ar‘c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:796:警告:未使用变量'NThist_ar’c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:797:警告:未使用变量'SThist_ar‘c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:798:警告:未使用变量'DThist_ar‘c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:817:警告:未使用变量'Ndiscont’c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:818:警告:未使用变量'Sdiscont‘c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:819:警告:未使用变量'Ddiscont’跟踪(最近调用)最后)：文件"C:\Documents and Settings\thao\Desktop\mackey-glass.py"，第33行，在文件"C:\Python27\lib\site-packages\pydelay_dde23.py"，第1120行，在运行编译器=‘gcc’文件第355行，在内联**kw)文件"C:\Python27\lib\site-packages\scipy\weave\inline_tools.py"，行482中，在compile_function verbose=verbose中，**kw)文件"C:\Python27\lib\site-packages\scipy\weave\ext_tools.py"，第367行，编译详细=详细，**kw)文件"C:\Python27\lib\site-packages\scipy\weave\build_tools.py"，第272行，在build_extension安装程序(name= module_name，ext_modules = ext，verbose=verb)文件"C:\Python27\lib\site-packages\numpy\distutils\core.py"，第186行中在安装程序中，返回old_setup(**new_attr)文件"C:\Python27\lib\distutils\core.py"，第169行，在安装程序中，"error：“+ str(msg) distutils.errors.CompileError: distutils.errors.CompileError:g++ -O2 -Wall -IC：\Python27 27\lib\site-packages\scipy\weave -IC:\Python27\include -IC:\Python27\PC -c -o c:\docume~1\thao\locals~1\temp\thong\python27_intermediate\compiler_a77d1132635f0379270bcb96a5e542fc\Release\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.o“失败，退出状态1在5.8s中完成，退出代码1

方法1(使用if else语句更新eqns，这是一个python )：

import numpy as np
import pylab as pl
from pydelay import dde23

eqn_1a='0.25 * x(t-tau) / (1.0 + pow(x(t-tau),10.0)) -0.1*x'
eqn_1b='0.45 * x'

eqns = { 'x' : eqn_1a if 'x>1' else eqn_1b}

dde = dde23(eqns=eqns, params={'tau': 15})
dde.set_sim_params(tfinal=1000, dtmax=1.0, AbsTol=10**-6, RelTol=10**-3)

histfunc = {'x': lambda t: 0.5 } 
dde.hist_from_funcs(histfunc, 51)
dde.run()

sol1 = dde.sample(2, 3, 0.1)
x1 = sol1['x']

方法2(向求解者提供c代码)：

eqns = { 'x' : 'f(x, tau)'}

# We can define a c function to be used in the equations
mycode = """
    double f(double x, double tau) {
        if (x>1){
            return (0.25 * x(t-tau) / (1.0 + pow(x(t-tau),10.0)) -0.1*x);
        }
        else{
            return (0.45 * x);
        }
     }
    """
dde = dde23(eqns=eqns, params={'tau': 15}, supportcode=mycode)

dde.set_sim_params(tfinal=1000, dtmax=1.0, AbsTol=10**-6, RelTol=10**-3)

histfunc = {'x': lambda t: 0.5 } 
dde.hist_from_funcs(histfunc, 51)
dde.run()

sol1 = dde.sample(1, 300, 1)
x1 = sol1['x']

Answer 1

第二种方法只需做一个小小的改动。延迟变量的值必须作为参数直接传递给函数。

eqns = { 'x' : 'f(x, x(t - tau))'}

mycode = """
double f(double x, double x_tau) {
    if (x > 1.0){
        return 0.25 * x_tau / (1.0 + pow(x_tau, 10.0)) -0.1*x;
    }
    else{
        return 0.45 * x;
    }
}
"""

虽然这是运行的，但请注意，它不能给出一个非常准确的解决方案。这是由于x=1的不连续性。更高级的求解器可以显式地处理这种不连续(参见http://www.radford.edu/~thompson/ffddes/)。

如果您想为了方便而坚持使用pydelay解决程序，或者只是为了获得一个概述，您可以尝试将最大步长dtmax设置得足够小，以减少错误。