标签: ode

y''' + 41y'' + 360y' + 900y = 600x' + 1200x;
y(0)= 2 ; y'(0)= 1 ; y''(0) = -0.05

如何使用ODE45函数求解该等式？

我试过这个:

==>
function dydt=f(t,y)

dydt = [y(2) ; y(3) ; -41*y(3)-360*y(2)- 900*y(1)]

==>
clear all;

timerange=[0 1.4]; %seconds
initialvalues=[2 1 -0.05];


[t,y]=ode45(@dydt, timerange, initialvalues)

plot(t,y(:,1));

但我需要将X部分放在等式中 - 我不知道如何......

我想从数字的角度理解混合建模(特别是状态事件)背后的一般思路(虽然我不是数学家:)).鉴于以下Modelica模型:

model BouncingBall
    constant Real g=9.81 
    Real h(start=1);
    Real v(start=0);
equation
    der(h)=v;
    der(v)=-g;
algorithm
    when h < 0 then
         reinit(v,-pre(v));
    end when;

end BouncingBall;

我理解的概念when和reinit.

1. when声明中的等式仅在条件成立时才有效吗？
2. 让我们假设球会完全击中地板2sec.由于我使用多步求解器,这意味着求解器"超过2秒",识别出h<0(假设在模拟time = 2.5sec中h = -0.7).这是什么意思"使用交叉函数搜索事件的时间？是否有简单的解释(示例)？
3. 解算器现在要回去了吗？采取较小的步长？
4. pre()在这种情况下,操作意味着什么？
5. noEvent():"表达式是字面意义而不是生成交叉函数.由于没有交叉函数,因此没有要求表达式可以在事件限制之外进行评估":这是什么意思？给出与弹跳球相同的例子:求解器在时间2.5检测到h = 0.7.有没有区别noEvent()？

我试图使用nlme和拟合一阶微分模型lsoda。这是基本思想：我首先定义允许生成微分方程解的函数：

library(deSolve)

ODE1 <- function(time, x, parms) {with(as.list(c(parms, x)), {
  import <- excfunc(time)
  dS <- import*k/tau - (S-yo)/tau 
  res <- c(dS)
  list(res)})}


solution_ODE1 = function(tau1,k1,yo1,excitation,time){
  excfunc <- approxfun(time, excitation, rule = 2)
  parms  <- c(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)
  xstart = c(S = yo1)
  out <-  lsoda(xstart, time, ODE1, parms)
  return(out[,2])
}

然后，根据两个ID的公式生成数据：

time <- 0:49
excitation <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))
simu_data <- data.frame(signal = c(solution_ODE1(3,2,0.1,excitation,time)+rnorm(length(time),0,0.1),
                                   solution_ODE1(3.2,1.5,0.3,excitation,time)+rnorm(length(time),0,0.1)),
                        time = rep(time,2),
                        excitation …

我有一个二阶微分方程，我想在 python 中解决它。问题是，对于其中一个变量，我没有初始条件，0而只有无穷大的值。谁能告诉我应该提供哪些参数scipy.integrate.odeint？能解决吗？

方程：

需要在时间方面找到 Theta。它的一阶导数在 处为零t=0。theta 未知，t=0但它在足够长的时间内变为零。其余的都是已知的。由于近似值I可以设置为零，因此消除了二阶导数，这将使问题更容易。

我正在尝试使用 deSolve 在 R 中解决一个简单的 ODE： dQ/dt = f(Q)*(P - E).整个过程是 Q 的时间序列。诀窍是 P 和 E 是数据本身的固定时间序列，因此 diff eq 仅在 Q 中有效。

用固定的时间步长迭代解决这个问题相对简单，但我试图找到一种在 R 中使用自适应时间步长的方法。因为 P 和 E 是离散的，连续的时间步长可能具有相同的 P 和 E 值，其中很好。我一直在玩 deSolve，但一直无法解决这个问题。理想情况下，我想使用标准的四阶 Runge-Kutta。

有任何想法吗？在 MATLAB 中做吗？

编辑可复制的例子。我希望能够使用可变时间步长 Runge-Kutta 4 方法进行此计算。我可以很容易地对固定时间步 rk4 进行编程，而不是自适应。

working <- structure(list(datetime = structure(c(1185915600, 1185919200, 
1185922800, 1185926400, 1185930000, 1185933600, 1185937200, 1185940800, 
1185944400, 1185948000, 1185951600), class = c("POSIXct", "POSIXt"
), tzone = "UTC"), p = c(0, 0, 0, 1.1, 0.5, 0.7, 0, 0, 1.3, 0, …

我正在尝试数值求解一个允许离散跳跃的 ODE。我正在使用 Euler 方法，并希望 Numba 的 jit 可以帮助我加快进程（现在脚本需要 300 秒才能运行，而我需要它运行 200 次）。

这是我简化的第一次尝试：

import numpy as np
from numba import jit

dt = 1e-5
T = 1
x0 = 1
noiter = int(T / dt)
res = np.zeros(noiter)

def fdot(x, t):
    return -x + t / (x + 1) ** 2

def solve_my_ODE(res, fdot, x0, T, dt):
    res[0] = x0
    noiter = int(T / dt)
    for i in range(noiter - 1):
        res[i + 1] = res[i] + dt * fdot(res[i], …

我正在寻找在满足特定条件时终止 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') …

我有一个包含 3 个边界条件的 3 个微分方程系统（从我相信的代码中可以很明显地看出）。我设法在 MATLAB 中用一个循环来解决它，以便在程序即将返回错误时一点一点地改变初始猜测而不终止程序。然而，在scipy's solve_bvp，我总是得到一些答案，虽然它是错误的。所以我一直在改变我的猜测（这一直在改变答案）并且我给出的数字与我从实际解决方案中得到的数字非常接近，但它仍然无法正常工作。代码是否存在其他问题，因此无法正常工作？我刚刚编辑了他们的文档代码。

import numpy as np
def fun(x, y):
    return np.vstack((3.769911184e12*np.exp(-19846/y[1])*(1-y[0]), 0.2056315191*(y[2]-y[1])+6.511664773e14*np.exp(-19846/y[1])*(1-y[0]), 1.696460033*(y[2]-y[1])))
def bc(ya, yb):
    return np.array([ya[0], ya[1]-673, yb[2]-200])
x = np.linspace(0, 1, 5)
#y = np.ones((3, x.size))
y = np.array([[1, 1, 1, 1, 1], [670, 670, 670, 670, 670], [670, 670, 670, 670, 670] ])
from scipy.integrate import solve_bvp
sol = solve_bvp(fun, bc, x, y)

实际解决方案如下图所示。

BVP 的 MATLAB 解决方案

我正在尝试做一个简单的谐振子示例，它将通过龙格-库塔四阶方法求解。待解的二阶常微分方程（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 = …

我有一个函数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 - …