标签: ode

我是Python初学者。我正在尝试切换 matlab 中的一些程序。我需要求解一个刚性颂方程，其输入都是矩阵。在matlab中我使用

[ttT,uT] = ode23s('SST',t,fT);

我喜欢求解涉及多个阈值的耦合微分方程组。通过查看 R 信息，我将 ODE 与根函数和事件函数结合使用。

通过各种示例，即温度模型，第 14 页http://cran.r-project.org/web/packages/diffEq/vignettes/ODEinR.pdf - 下面粘贴的代码 - 我能够让我的模型发挥作用在一个阈值上，即找到一个变量的根/达到阈值会触发一个事件。

library(deSolve)\nyini <- c(temp = 18, heating_on=1)\n\ntemp <- function(t,y, parms) {\n  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)\n  dy2 <- 0\n list(c(dy1, dy2))\n}\n\nrootfunc <- function(t,y,parms) c(y[1]-18, y[1]-20)\n\neventfunc <- function(t,y,parms) {\n  y[1] <- y[1]\n  y[2] <- ! y[2]  \n return(y)\n}\n\ntimes <- seq(from=0, to=20, by=0.1)\nout <- lsode(times=times, y=yini, func = temp, parms = NULL, \n         rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))\nplot(out, lwd=2)\nattributes(out)$troot\n

该示例还表明，不同的根可以触发相同的事件函数（y[1] \xe2\x80\x93 …

我正在寻找在满足特定条件时终止 MATLAB ode 的方法。我在本主题MatLab ODE 启动/停止条件中找到了答案 ，其中讨论了“事件”的使用。然而，这适用于 ode45，当我尝试将“事件”与 ode15i 一起使用时，它根本不起作用，并且 MATLAB 显示错误。

我试图通过简单的例子来学习这一点，并求解一个简单的微分方程组，如下所示。

dx/dt = 5x + 3y；dy/dt = x + 7y；我使用 ode45 解决了它们，并尝试使用 ode15i 执行相同的操作，但它不起作用。下面给出的是我的代码。

与 ode45

    function start_stop_test_ode45
    y0 = [5;1];
    tv = linspace(0,2,100);
    options = odeset('Events',@events);    
    f = @(t,y) [5*y(1) + 3*y(2);y(1) + 7*y(2)];
    [t,Y] = ode45(f,tv,y0,options);
    xNI = Y(:,1);
    yNI = Y(:,2);
    xCF   =  3*exp(4*t) + 2*exp(8*t);
    yCF   = -1*exp(4*t) + 2*exp(8*t);
    % Here we plot all the graphs
    figure(1)
    plot(t,xNI,'--k',t,xCF,'r','Linewidth',1.75)
    xlabel('t (s)')
    ylabel('x') …

这里scipy.integrate.odeint用六个不同的标准颂歌问题来调用rtol= atolfrom 1E-06to 1E-13。我查看了所有较大公差的结果之间的最大差异减去最小公差的结果，以获得某种“错误”的表示。我很好奇为什么对于给定的容差，一个问题 (D5) 给出的错误比另一个问题 (C1) 严重一百万倍，即使步数范围相当窄（在 10 倍之内）。

脚本中给出了颂歌问题的引用。所有问题都相当正常化，所以我正在进行rtol类似的治疗atol。

重申一下 - 我的问题是，为什么不同问题之间的误差几乎相差一个因子1E+06，尽管误差随容差而变化。当然，C1 是“最柔和的”，D5 在“近日点”处具有戏剧性的峰值，但我认为例程会在内部调整步长，以便误差相似。

编辑：我添加了“错误”的时间演变，这可能会带来一些启发。

# FROM: "Comparing Numerical Methods for Ordinary Differential Equations"
# T.E. Hull, W.H. Enright, B.M. Fellen and A.E. Sedgwidh
# SIAM J. Numer. Anal. vol 9, no 4, December 1972, pp: 603-637

def deriv_B1(y, x):
    return [2.*(y[0]-y[0]*y[1]), -(y[1]-y[0]*y[1])] # "growth of two conflicting populations"

def deriv_B4(y, x):
    A = …

我正在开发使用 scipy 的 complex_ode 集成 ODE 的代码，其中被积函数包括傅立叶变换和作用于大量复数值的指数运算符。

为了优化性能，我已经对此进行了分析，发现主要瓶颈是（在使用 PyFFTW 等优化 FFT 之后）：

val = np.exp(float_value * arr)

我目前正在使用 numpy，我理解它称为 C 代码 - 因此应该很快。但是有什么办法可以进一步提高性能吗？

我已经研究过使用 Numba，但由于我的主循环也包含 FFT，我认为它无法编译（nopython=True 标志会导致错误），因此，我怀疑它没有任何收益。

这是我要优化的代码的测试示例：

arr = np.random.rand(2**14) + 1j *np.random.rand(2**14)
float_value = 0.5
%timeit np.exp(float_value * arr)

任何建议欢迎感谢。

我有用于绘制以下捕食者猎物模型的代码：

dx/dt = x ? xy, dy/dt = ?y + xy

from pylab import *
xvalues, yvalues = meshgrid(arange(0, 3, 0.1), arange(0, 3, 0.1))
xdot = xvalues - xvalues * yvalues
ydot = - yvalues + xvalues * yvalues
streamplot(xvalues, yvalues, xdot, ydot)
show()

但我不确定如何使用这些函数来绘制相平面（使用流图）来模拟钟摆运动，定义为

d^2?/dt^2 = (?g/L)sin(?)

如何使用 matplotlib 和 pylab 实现此模型以生成相平面？

我正在尝试做一个简单的谐振子示例，它将通过龙格-库塔四阶方法求解。待解的二阶常微分方程（ODE）及初始条件为：

y'' + y = 0

y(0) = 0 且 y'(0) = 1/pi

范围在 0 到 1 之间，有 100 级。我将二阶 ODE 分成两个一阶 ODE，使用 u 作为辅助变量：

y' = u

你' = -y

解析解为正弦 y(x) = (1/pi)^2 sin(pi*x)。

我的Python代码如下：

from math import pi
from numpy import arange
from matplotlib.pyplot import plot, show

# y' = u
# u' = -y

def F(y, u, x):
    return -y

a = 0
b = 1.0
N =100
h = (b-a)/N

xpoints = arange(a,b,h)
ypoints = …

我正在尝试编写一个函数，该函数将使用 4 阶隐式 Runge-Kutta 方法 (IRK) 求解 ODES 系统，但我无法正确定义循环。这里我们定义 IRK

任何建议将不胜感激！

function [tout,yout] = IRK4Solver(f,t,y0) 
t = t(:); % ensure that t is a column vector
N = length(t); 
h = (t(end)-t(1))/(N-1); % calculate h by assuming t gridpoints are equidistant
d = length(y0); % defines the size of y0
y0 = y0(:); % ensures that y0 is a column vector
y = zeros(d,N);% allocate the output array y
y(:,1) = y0; % assign y0 to the first column of …

我是 Julia 编程的新手，我设法解决了一些一阶 ODE，但是当我想转向二阶时，我不知道如何使用求解器来实现所需的方程。

我想解这个方程

y" + y = 0

具有初始条件

y(0) = 3
y'(0) = -0.5

我怎样才能做到这一点？

我有一个函数d\xcf\x86/dt = \xce\xb3 - F(\xcf\x86)（其中F(\xcf\x86)-- a 是2\xcf\x80周期函数）和该函数的图形F(\xcf\x86)。

我需要创建一个程序，输出 6 个\xcf\x86(t)不同值\xce\xb3( \xce\xb3 = 0.1、0.5、0.95、1.05、2、5)和 的图t\xe2\x88\x88[0,100]。

这是函数的定义F(\xcf\x86)：

      -\xcf\x86/a - \xcf\x80/a,    if \xcf\x86 \xe2\x88\x88 [-\xcf\x80, -\xcf\x80 + a]\n      -1,            if \xcf\x86 \xe2\x88\x88 [-\xcf\x80 + a, - a] \nF(\xcf\x86) = \xcf\x86/a,          if \xcf\x86 \xe2\x88\x88 [- a, a]\n       1,            if \xcf\x86 \xe2\x88\x88 [a, \xcf\x80 - …