标签: ode

我有以下 Matlab ODE 代码：

[t,y,~,~,ie] = ode23tb(@(t,y) RHSODE(t,y),[0,t_end], [i0;v0],options);

我希望 ODE 求解器还可以给出结果 z，它是 y 和 dy/dt 的函数，这样 z = f(y,dy/dt)。

有谁知道如何将这样的 z 添加到求解器的输出中？

ode我正在尝试使用包中的函数求解一阶微分方程deSolve。问题如下：药物在某些时间（输注时间）以恒定的输注速率给药，并以一阶速率消除。因此，该过程可以描述为：

if(t %in% Infusion_times){Infusion <- Infusion_rate} else{Infusion <- 0}
    dC <- -Ke*C + Infusion

其中t是时间，Infusion_times是包含给药时间的向量，C是药物的量，Ke是其消除常数，Infusion是在输注时取输注速率值的变量，否则值为 0。因此，假设我们想要注射 3 剂，从时间 0、24 和 40 开始，每次输注持续两个小时，假设我们想要deSolve每 0.02 个时间单位计算一次答案。例如，我们想要deSolve求解 0 到 48 之间时间的微分方程，步长为 0.02 倍单位。所以我们指定函数时间的向量ode将是：

times <- seq(from = 0, to = 48, by = 0.02)

输注时间由下式给出：

Infusion_times <- c(seq(from = 0, to = 2, by = 0.02), seq(from = 24, to = 26, by = …

我想用数值方法求解以下常微分方程：

其中 W 是白噪声，delta 是狄拉克 delta 函数。如果系统只有 p 和 v 中的项，那么用 R 编写代码会很容易。事实上，这是我的工作：

# parameters:
parameters <- c(alpha = 1,
                beta = 2,
                kappa2 = 1,
                kappa3 =-1)

# initial conditions:
state <- c(v=0.8,p=0.5)

Lorenz<-function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    # rate of change
      dv <- -2*alpha*v-beta^2*p
      dp <- v
      # return the rate of change
      list(c(dv,dp))
  }) # end with(as.list ...
}


times <- seq(0, 100, by = 0.01)


library(deSolve)
out <- ode(y = state, times = times, func …

我正在尝试复制此处发布的湖泊食物网络模型。该模型代表两条食物链（沿海与远洋），由顶级捕食者（鱼类）连接。我已经对模型进行了编码，但是当我在 2-3 个时间步长后运行它时，模型会生成NaN. 我已经多次检查我的代码，寻找括号等问题，但找不到问题。

如果我将fish初始丰度设置为 0，模型就会运行，所以我认为问题一定出在模型的鱼组件上。

以下是方程式：

Ap = 中上层资源，Z = 中上层浮游动物，Pp = 中上层捕食者，F = 鱼类，Al = 沿岸资源，I = 无脊椎动物，Pl = 沿岸捕食者。

这是我对模型进行编码的尝试：

library(deSolve)

# define the model
vade_2005_model <- function(Time, State, Pars){
  
  with(as.list(c(State, Pars)), {

# pelagic components -----------------------------------------------

# resource
pel_res_dt <- (rPel * AP * (1 - (AP/(KT * q)))) - (aZP * ((Z * AP)/(AP + bZP)))

# zooplankton
pel_zoo_dt <-  (ef * aZP * ((Z * AP)/(AP + bZP))) …

让我们说我有一个简单的逻辑方程

dx/dt = 2ax(1-x/N) 其中N是承载能力,a是一些增长率,a和N都是我想要变化的参数.

所以我想要做的是绘制我的固定点和两个参数的3D图形.

我理解如何找到单个参数的固定点.

这是我的示例代码

function xprime = MyLogisticFunction(t,X) %% The ODE

% Parameters
N = 10 % Carrying Capacity
a = 0.5 % Growth Rate

x1prime = 2*a*X(1)*(1 - X(1)/N );

xprime = [x1prime ]';

end

接下来是我的解算器

% Initial Number 
x0 = 0.4;

%Time Window
tspan=[0 100];

[t,x]=ode45(@MyLogisticFunction,tspan,x0);

clf

x(end,1) % This gives me the fixed point for the parameters above.

所以我真正的问题是,如何在两个函数之间放置一个for循环,这允许我改变a和N,这样我就可以绘制出 …

我想解决：

[\mathbf{M} \ddot{ \mathbf{U} }+ \mathbf{C} \dot{ \mathbf{U} }+ \mathbf{K} \mathbf{U} = \mathbf{P}(t) ]

或者，以状态空间形式：

[\dot{\mathbf{Y}}=f(\mathbf{Y},t)]

在哪里：

[\mathbf{Y} = \left[ \begin{array}{ c} \mathbf{U} \ \dot{ \mathbf{U} \end{array} \right] ]

和：

[f( \mathbf{Y} ,t)= \left[ \begin{array}{ c} \dot{ \mathbf{U} }\ \mathbf{M}^{-1} \mathbf{P} (t )- \mathbf{M} ^{-1} \mathbf{C} \dot{ \mathbf{U} }- \mathbf{M} ^{-1} \mathbf{K} \mathbf{U} \end{数组} \right] ]

我在 Julia 中尝试了以下代码，使用

\mathbf{M} = \left[ \begin{array}{ cc} 2&0\ 0&1 \end{array} \right];

\mathbf{C} = \left[ \begin{array}{ cc} 0&0\ 0& 0 \end{array} \right];

\mathbf{K} = …

在 Julia 中，我想解决具有外部强迫的 ODE 系统，g1(t), g2(t)例如

dx1(t) / dt = f1(x1, t) + g1(t)
dx2(t) / dt = f2(x1, x2, t) + g2(t)

与从文件中读入的强制。

我正在使用这项研究来学习 Julia 和微分方程包，但我很难找到正确的方法。

我可以想象使用 acallback可以工作，但这似乎很麻烦。

您知道如何实施这种外部强迫吗？

在 Sympy 中 ODE 的定义（和解决方案）中是否有充分的理由使用Derivative而不是？似乎做得很好：diffdiff

我正在研究非线性微分方程.我所做的是平均超过100个不同初始条件值的位置.

我在gsl中使用了odeiv.对于每个初始值,时间范围是4*10 ^ 7.然而,程序杀死,一旦我设置了10个不同的初始条件和10 ^ 6的时间范围.这是一种限制.

我的电脑有8个内核和16GB内存.我认为这不是那么大.

我将把编码的一部分.有人帮我这个吗？谢谢.

long long int i, j, k;
double const y_i = 0, dydt_i = 0;
double n = 10;
long long int tmax = 1000000;
double A[tmax];

for (j=0; j < n; j++)
{
    double y_j = y_i + 0.001*j;
    double dydt_j = dydt_i;
    t = 0.0;
    double y[2] = {y_j, dydt_j};
    gsl_odeiv2_system sys = {func, jac, 2, &params};
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6,     0.0);

  for (i=0; i< tmax; …

我想解决以下颂歌

K T + C T' = Q

给定的示例数据是我下面的代码

import numpy as np
import scipy as sp
# Solve the following ODE
# K*T + C*T' = Q
# T' = C^-1 ( Q - K * T )
T_start=sp.array([ 151.26,  132.18,  131.64,  146.55,  147.87,  137.87])

K = sp.array([[-0.01761969,  0.02704873,  0.00572222,  0.        ,  0.        ,
         0.        ],
       [ 0.02704873, -0.03546941,  0.        ,  0.        ,  0.00513177,
         0.        ],
       [ 0.00572222,  0.        ,  0.03001858, -0.04752982,  0.        ,
         0.02030505],
       [ 0.        ,  0. …

我希望有人可以帮助我编写一些简单的代码来压缩向量。

我目前正在研究 ODE 疾病模型。我有一个代表受感染马匹数量的向量。马匹在离开车厢之前必须经过“x”个阶段。该模型是在网格系统上建立的，因此对于每次模拟，我都有“n”个在向量中表示的单元格。所以向量长度是x*n。

我想总结每个单元格所有阶段的马匹数量，因此结果是长度为“n”的向量。那么例如...

n <- 5
x <- 3

inf <- c(1,2,3,
         4,5,6,
         7,8,9,
         10,11,12,
         13,14,15)

outcome <- c(6,
             15,
             24,
             33,
             42)

理想情况下，我希望不必将其放入某种循环中，因为我的代码计算量已经非常大，而且我无法增加更多的运行时间。