我有两个 numpy 数组：9x9 和 9x1。我想求解离散时间点的微分方程，但在让 ODEInt 工作时遇到问题。我确实不确定我是否在做正确的事。
对于 Mathematica，方程是：
Solution = {A[t]} /. NDSolve[{A'[t] == Ab.A[t] && A[0] == A0}, {A[t]}, {t, 0, .5}, MaxSteps -> \[Infinity]];
time = 0.25;
increment = 0.05;
MA = Table[Solution, {t, 0, time, increment}];
其中 Ab 是 9x9 矩阵，A0 是 9x1 矩阵（初始）。在这里，我解决了时间问题，生活很美好。
在Python实现中，我有以下代码，它给了我错误的答案：
from scipy.integrate import odeint
from numpy import array, dot, pi
def deriv(A, t, Ab):
return dot(Ab, A)
def MatrixBM3(k12,k21,k13,k31,k23,k32,delta1,delta2,delta3,
w1, R1, R2):
K = array([[-k21 -k23, k12, k32, 0., 0., 0., …SciPy 可以通过 scipy.integrate.odeint 或其他包求解 ode 方程，但它在函数完全求解后给出结果。但是，如果 ode 函数非常复杂，程序将花费大量时间（一两天）才能给出整个结果。那么我如何监控求解方程的步骤（当方程尚未完全求解时打印出结果）？
假设我有以下boost::odeint代码:
#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
using namespace std;
using namespace boost::numeric::odeint;
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
typedef boost::array< double , 3 > state_type;
void lorenz( const state_type &x , state_type &dxdt , double t ){
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = -b * x[2] + x[0] * x[1]; …我在 Boost 中使用 odeint 库。使用integrate_adaptive函数时，结果符合预期。但是，使用Integrate_times 时，ODE 会在非常不同的时间进行计算，这些时间超出了积分范围。这对我来说是一个问题，因为我的 ODE 没有针对某些正在评估的值进行定义。
下面的代码演示了这个问题。计算 ODE 的 x 值将打印到屏幕上。
#include <iostream>
#include <complex>
#include <vector>
#include <boost/numeric/odeint.hpp>
struct observe
{
std::vector<std::vector<std::complex<double> > > & y;
std::vector<double>& x_ode;
observe(std::vector<std::vector<std::complex<double> > > &p_y, std::vector<double> &p_x_ode) : y(p_y), x_ode(p_x_ode) { };
void operator()(const std::vector<std::complex<double> > &y_temp, double x_temp)
{
y.push_back(y_temp);
x_ode.push_back(x_temp);
}
};
class Direct
{
std::complex<double> alpha;
std::complex<double> beta;
std::complex<double> R;
std::vector<std::vector<std::complex<double> > > H0_create(const double y);
public:
Direct(std::complex<double> p_alpha, std::complex<double> p_beta, double p_R) : alpha(p_alpha), beta(p_beta), …我想runge_kutta4在odeint C++库中使用方法.我在Matlab中解决了这个问题.在Matlab我以下代码来解决x'' = -x - g*x',具有初始值的x1 = 1,x2 = 0是如下
main.m
clear all
clc
t = 0:0.1:10;
x0 = [1; 0];
[t, x] = ode45('ODESolver', t, x0);
plot(t, x(:,1));
title('Position');
xlabel('time (sec)');
ylabel('x(t)');
ODESolver.m
function dx = ODESolver(t, x)
dx = zeros(2,1);
g = 0.15;
dx(1) = x(2);
dx(2) = -x(1) - g*x(2);
end
我已经安装了odeint库.我的使用代码runge_kutta4如下
#include <iostream>
#include <boost/numeric/odeint.hpp>
using namespace std;
using namespace boost::numeric::odeint;
/* The type …我正在尝试使用numba来提升scipy.integrate.odeint的python性能.为此,我必须使用@ nb.jit(nopython = True)来定义ODE系统的函数.但是,这个函数必须在我的程序中将另一个python类实例作为参数.我还必须使用适当的规格与@ nb.jitclass(spec)一起使用jit.这很好,直到我发现一个严重的问题,当类的规范包括另一种类的实例作为其方法.我的代码如下.
import numba as nb
from scipy.integrate import odeint
spec1=[("hi", nb.i4)]
@nb.jitclass(spec1)
class Hi(object):
def __init__(self):
self.hi = 0
spec2=[("dummy", nb.i4), ("dummy1", nb.i4)]
@nb.jitclass(spec2)
class Dummy(object):
def __init__(self, anotherClassInstance):
self.dummy = 0
self.dummy1 = anotherClassInstance
class A:
def __init__(self, someClassInstance):
self.a1 = someClassInstance
def odeSystem(self, x, t):
return _odeSystem(x, t, self.a1)
def odeSolve(self, iValues, ts):
sol = odeint(self.odeSystem, iValues, ts)
return sol
@nb.jit(nopython=True)
def _odeSystem(x, t, someClassInstance):
return 1-x
if __name__ == "__main__":
c = Hi() …我尝试在Debian Squeeze g ++ 4.4上的boost_1_54_0中运行odeint示例
Lorenz系统工作正常,但简单的1d ode:
#include <iostream>
#include <boost/numeric/odeint.hpp>
using namespace std;
using namespace boost::numeric::odeint;
/* we solve the simple ODE x' = 3/(2t^2) + x/(2t)
* with initial condition x(1) = 0.
* Analytic solution is x(t) = sqrt(t) - 1/t
*/
void rhs( const double x , double &dxdt , const double t )
{
dxdt = 3.0/(2.0*t*t) + x/(2.0*t);
}
void write_cout( const double &x , const double t …我需要解决具有复杂初始值的复杂域定义的ODE系统.scipy.integrate.odeint不适用于复杂系统.我强调在实部和虚部切割我的系统并分别求解,但我的ODE系统的rhs涉及因变量本身和它们的复共轭之间的产物.哇,我这样做吗？这是我的代码,我尝试在Re和Im部分中破坏RHS,但我不认为解决方案是否相同,因为复杂数字之间的内部产品我不会打破它.在我的脚本中,u1是一个(非常)长的复杂函数,比如u1(Lm)= f_real(Lm)+ 1j*f_imag(Lm).
from numpy import *
from scipy import integrate
def cj(z): return z.conjugate()
def dydt(y, t=0):
# Notation
# Dependent Variables
theta1 = y[0]
theta3 = y[1]
Lm = y[2]
u11 = u1(Lm)
u13 = u1(3*Lm)
zeta1 = -2*E*u11*theta1
zeta3 = -2*E*3*u13*theta3
# Coefficients
A0 = theta1*cj(zeta1) + 3*theta3*cj(zeta3)
A2 = -zeta1*theta1 + 3*cj(zeta1)*theta3 + zeta3*cj(theta1)
A4 = -theta1*zeta3 - 3*zeta1*theta3
A6 = -3*theta3*zeta3
A = - (A2/2 + A4/4 + A6/6)
# RHS vector components
dy1dt = …我对编码完全不熟悉,我想用数字方法解决这5个微分方程.我拿了一个python模板并将它应用到我的案例中.这是我写的简化版:
import numpy as np
from math import *
from matplotlib import rc, font_manager
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#Constants and parameters
alpha=1/137.
k=1.e-9
T=40.
V= 6.e-6
r = 6.9673e12
u = 1.51856e7
#defining dy/dt's
def f(y, t):
A = y[0]
B = y[1]
C = y[2]
D = y[3]
E = y[4]
# the model equations
f0 = 1.519e21*(-2*k/T*(k - (alpha/pi)*(B+V))*A)
f1 = (3*B**2 + 3*C**2 + 6*B*C + 2*pi**2*B*T + pi**2*T**2)**-1*(-f0*alpha/(3*pi**3) - 2*r*(B**3 …我想知道，求解微分方程之间ODEINT和solve_ivp求解微分方程有什么区别。它们之间可能有什么优点和缺点？
f1 = solve_ivp(f, [0,1], y0) #y0 is the initial point
f2 = odeint(f, y0, [0, 1], args=(a, b)) # a and b are arguments of function f
谢谢
odeint在Numpy 模拟过程中保存中间变量的最简单方法是什么？
例如:
def dy(y,t)
x = np.rand(3,1)
return y + x.sum()
sim = odeint(dy,0,np.arange(0,1,0.1))
x在模拟过程中保存存储数据的最简单方法是什么？理想情况下,在t传递给的参数中指定的点odeint.
我对C++很新,我试图通过Rcpp来加速我的R代码.
下面的代码从t0到t1集成 - 这是在"lorenz"函数中完成的.Test4使用"lorenz""计数"次数进行整合.但是,在时间"t1",系统的状态在系统重新运行之前在"write_lorenz"中被修改,这就是问题所在.如果我通过从R调用test4反复运行相同的程序,打印到屏幕总是产生相同的结果,但是,我返回的矩阵"u"没有,并且似乎最终收敛到任何"t1"是什么问题.
我的输入值没有改变所以我想知道是否有内存泄漏,或者是否有其他事情发生,如何解决它.另外我想知道我的"u"初始化是否不正确我应该使用"new"命令.
我试过的是NumericMatrix*u = NULL;*u =新的NumericMatrix; 然后我尝试访问矩阵的元素为*u(1,2),但是以这种方式访问元素会导致错误,说你不是函数.
任何帮助将不胜感激
我从以下网站修改了此代码
http://headmyshoulder.github.io/odeint-v2/examples.html
所以我可以和Rcpp一起使用它
//###################################R Code ###############################
library(Rcpp)
sourceCpp("test4.cpp")
sigma <- 10.0
R <-28.0
b <- 8.0 / 3.0
a2 <- c(10.0 , 1.0 , 1.0) #initial conditions X0,I0,Y0
L2 <- c(0.0 , 2.0 , 0.1) #initial time, kick time, error
counts <- 2
kick <-1.0; # kick size
pars <-c(sigma,R,b,kick)
test4(a=a,L2=L2,counts=counts,pars= pars)
// C ++ code
//[[Rcpp::depends(BH)]]
//[[Rcpp::depends(RcppEigen)]]
//[[Rcpp::plugins("cpp11")]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <math.h>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp> …你能为我提供执行与数字融合的一个简单的例子odeint的C++？
我想使用方便的集成功能,记录为:
integrate( system , x0 , t0 , t1 , dt )
另外我不确定,如果可能的话,如何传递它而不是函数或仿函数,类方法.