标签: ode

我是 Octave 的新手，所以在转向更复杂的项目之前，我试图让一些简单的示例发挥作用。

我正在尝试解决 ODE 问题dy/dx = a*x+b，但没有成功。这是代码：

%Funzione retta y = a*x + b. Ingressi: vettore valori t; coefficienti a,b
clear all;
%Inizializza argomenti
b = 1;
a = 1;
x = ones(1,20);
function y = retta(a, x, b) %Definisce funzione
y = ones(1,20);
y = a .* x .+ b;
endfunction
%Calcola retta
x = [-10:10];
a = 2;
b = 2;
r = retta(a, x, b)
c = b;
p1 = (a/2)*x.^2+b.*x+c  %Sol. analitica di …

考虑到摆锤初始角度(x),重力加速度(g),线长度(l)和时间步长(h),我试图求解摆动运动的微分方程.我用Euler方法试过这个,一切都没问题.但现在我要使用在GSL中实现的Runge-Kutta方法.我试图从gsl手册中学习它,但我遇到了一个问题.钟摆不想停下来.假设我以初始角度1 rad开始它,无论它已经做了多少波动,它始终具有1弧度的峰值倾斜.这是我用来给GSL的公式和函数:

x''(t) + g/l*sin(x(t)) = 0

改造它:

x''(t) = -g/l*sin(x(t))

和分解:

y(t) = x'(t)
y'(t) = -g/l*sin(x(t))

这是代码片段,如果这还不够,我可以发布整个程序(它不会太长),但也许这里的问题在某处:

    int func (double t, const double x[], double dxdt[], void *params){
        double l = *(double*) params;
        double g = *(double*) (params+sizeof(double));
        dxdt[0] = x[1];
        dxdt[1] = -g/l*sin(x[0]);
        return GSL_SUCCESS;
    }

参数g并l正确传递给函数,我已经检查过了.

我想将n阶常微分方程简化为一阶方程组.这是为数值分析做准备.我使用Sympy和Sagemath作为计算机代数,但我没有找到任何函数来进行这种类型的缩减.我不确定是否有其他人可以指出Sympy或Sagemath中是否存在此功能.

这方面的一个例子是减少等式:

x''' - 2x'' + x' = 0 

到一阶方程组:

[0  1 0 
 0  0 1 
 0 -1 2]

我有以下Lotka-Volterra模型

dN1/dt = N1(1-N1-0.7N2)

dN2/dt = N2(1-N2-0.3N1)

N旁边的1和2是下标.

我想用SciPy来解决这个问题并将结果可视化.我想用y轴上的N2和N1上的N1绘制一个图.如果在第一个等式中将N1设置为零,则得到N2 = 1/0.7,如果在第二个等式中将N2设置为零,则得到N1 = 0.3/1.这两条线假设相交.我如何在Python中执行此操作？

我在线阅读本教程(幻灯片6到16).这就是我到目前为止所拥有的.

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def derivN1(y,t):
    yprime=np.array([1-0.7y[0]])
    return yprime

def derivN2(y,t):
    yprime=np.array([1-0.3y[0]])
    return yprime

start=0
end=1
numsteps=1000
time=np.linspace(start,end,numsteps)
y0=np.array([10])

yN1=integrate.odeint(derivN1,y0,time)
yN2=integrate.odeint(derivN2,y0,time)

plt.plot(time,yN1[:])
plt.plot(time,yN2[:])

但情节不正确.更新:我认为我使用了错误的方法.我正在阅读另一个在线教程.我会更多地解决这个问题.在此期间,如果有人知道如何解决它,请告诉我.

对于二阶ODE（Python中的dopri5方法），以下代码始终会导致错误：C:\Users\MY\Anaconda3\lib\site-packages\scipy\integrate\_ode.py:1019: UserWarning: dopri5: larger nmax is needed self.messages.get(idid, 'Unexpected idid=%s' % idid))。我更改了参数，但似乎无济于事。即使设置nsteps=100000也不起作用。还有其他方法可以解决这个问题，而不仅仅是增加nsteps？

from scipy.integrate import ode
import numpy as np

def fun(t, y):
    return np.array([y[1], -3/t*y[1] + 7/(t**6)*y[0]])

yinit = np.array([0.01, 0.2])

dt = 0.01
t_stop = 2

solver = ode(fun).set_integrator('dopri5', nsteps=100000).set_initial_value(yinit)
solver.t = 0.001
t_RK4_sci = [0]
x_RK4_sci = [yinit]
while solver.successful() and solver.t < t_stop:
    solver.integrate(solver.t+dt, step=True)
    t_RK4_sci.append(solver.t)
    x_RK4_sci.append(solver.y)
t_RK4_sci = np.array(t_RK4_sci)
x_RK4_sci = np.array(x_RK4_sci)

我是python的新手。我有一个简单的微分系统，它由两个变量和两个微分方程和初始条件组成x0=1, y0=2：

dx/dt=6*y
dy/dt=(2t-3x)/4y

现在我正在尝试解决这两个微分方程，我选择odeint. 这是我的代码：

import matplotlib.pyplot as pl
import numpy as np
from scipy.integrate import odeint

def func(z,b):
    x, y=z
    return [6*y, (b-3*x)/(4*y)]    

z0=[1,2]
t = np.linspace(0,10,11)
b=2*t
xx=odeint(func, z0, b)
pl.figure(1)
pl.plot(t, xx[:,0])
pl.legend()
pl.show()

但结果不正确，并出现错误信息：

Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.

我不知道我的代码有什么问题，我该如何解决。任何帮助都会对我有用。

在发布实际代码之前,让我告诉您我的计算机的处理器和内存信息是好的:

昨天我发布了关于Lorenz方程(混沌理论中的经典方程),其中一个伟大的人帮助我并展示解决方案,这里是:

function f=lorenz(t,x,a,b,c)
        % solve differential equation like this
        %dx/dt=a*(y-x)
        %dy/dt=-x*z+b*x-y
        %dz/dt=xy-c*z/3
        f=zeros(3,1);% preallocate result
        f(1)=a*(x(2)-x(1));
        f(2)=-x(1)*x(3)+b*x(1)-x(2);
        f(3)=x(1)*x(2)-c*x(3)/3;
        end

和测试程序(脚本):

% test program
x0=[-2 -3.5 21];% initial  point
a=input(' enter first coefficient : ');
b=input(' enter second coefficient: ');
c=input(' enter third coefficient : ');
[t,x] = ode45(@(t,x) lorenz(t,x,a,b,c),[0 10],x0);
plot(t,x(:,1),'r'); 
title(' solution of x part');
grid on

但在运行这些线后

test_program
 enter first coefficient : 10
 enter second coefficient: 28
 enter third coefficient : -8

它还在运行,他说在他的个人电脑上需要2秒钟,所以发生了什么事真的很奇怪？为什么不在我的电脑上编译？即使你看到它我的笔记本电脑有很好的参数,请帮助我 - 即使现在它正在运行所以我应该取消使用ctrl-c.

我很清楚如何以scipy.integrate.ode.integrate(t)最简单的形式使用该函数,但API读取它还需要两个可选参数,即step和relax.当前文档没有关于这些参数的信息,也没有在示例中使用它们.我想知道,他们做了什么以及他们有用的情况是什么？

我的 ode 求解器遇到了一些问题，我正在尝试解决 SEIR 问题，尽管我的代码基于非常相似的代码，但我不断收到相同的错误。我的代码是：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Total population, N.
N1 = 55600
# Initial number of infected and recovered individuals, I0 and R0.
I10, R10, E10 = 1, 0, 0
# Everyone else, S0, is susceptible to infection initially.
S10 = N1 - I10 - R10 - E10
# parameters
B = 0.05
a = 0.000001
d = 0.0167
g = 0.0167
z = 0.0167
M = 100000 …

我尝试在 Julia 文档中使用这个示例。我的尝试是将细胞分成两部分，每部分的蛋白质含量都是一半，所以我设置了 Theta=0.5。然而，情节是这样的：

很明显，每次达到目标蛋白质量时，细胞数量都会增加一倍，同时，因为它们是相等的。我怎么能画这个？我也不明白为什么在下面的情况下单元格的数量停在 3。