Решение динамической системы ODE с помощью Python - PullRequest
Новая
       54

Решение динамической системы ODE с помощью Python

0 голосов
/ 18 ноября 2018

Я пытаюсь решить динамическую систему с тремя переменными состояния V1, V2, I3 и затем построить их в 3D-графике. Мой код пока выглядит следующим образом: 

from scipy.integrate import ode
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import math
def ID(V,a,b):
    return a*(math.exp(b*V)-math.exp(-b*V))


def dynamical_system(t,z,C1,C2,L,R1,R2,R3,RN,a,b):

    V1,V2,I3 = z
    f = [(1/C1)*(V1*(1/RN-1/R1)-ID(V1-V2,a,b)-(V1-V2)/R2),(1/C2)*(ID(V1-V2,a,b)+(V1-V2)/R2-I3),(1/L)*(-I3*R3+V2)]
    return f 

# Create an `ode` instance to solve the system of differential
# equations defined by `dynamical_system`, and set the solver method to 'dopri5'.
solver = ode(dynamical_system)
solver.set_integrator('dopri5')

# Set the initial value z(0) = z0.
C1=10
C2=100
L=0.32
R1=22
R2=14.5
R3=100
RN=6.9
a=2.295*10**(-5)
b=3.0038
solver.set_f_params(C1,C2,L,R1,R2,R3,RN,a,b)
t0 = 0.0
z0 = [-3, 0.5, 0.25] #here you can set the inital values V1,V2,I3
solver.set_initial_value(z0, t0)


# Create the array `t` of time values at which to compute
# the solution, and create an array to hold the solution.
# Put the initial value in the solution array.
t1 = 25
N = 200 #number of iterations
t = np.linspace(t0, t1, N)
sol = np.empty((N, 3))
sol[0] = z0

# Repeatedly call the `integrate` method to advance the
# solution to time t[k], and save the solution in sol[k].
k = 1
while solver.successful() and solver.t < t1:
    solver.integrate(t[k])
    sol[k] = solver.y
    k += 1


xlim = (-4,1)
ylim= (-1,1)
zlim=(-1,1)
fig=plt.figure()
ax=fig.gca(projection='3d')
#ax.view_init(35,-28)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_zlim(zlim)

print sol[:,0]
print sol[:,1]
print sol[:,2]
ax.plot3D(sol[:,0], sol[:,1], sol[:,2], 'gray')

plt.show()

Печать массивов, которые должны содержать решения sol [:, 0] и т. Д., Показывает, что, по-видимому, он постоянно заполняет его начальным значением. Кто-нибудь может помочь? Спасибо!

1 Ответ

0 голосов
/ 18 ноября 2018

Использование from __future__ import division.

Я не могу воспроизвести вашу проблему: я вижу постепенное изменение с -3 до -2.46838127, с 0,5 до 0,38022886 и с 0,25 до 0,00380239 (с резкимизменить с 0,25 до 0,00498674 на первом этапе).Это с Python 3.7.0, NumPy версии 1.15.3 и SciPy версии 1.1.0.

Учитывая, что вы используете Python 2.7, здесь может быть целочисленное деление.Довольно много ваших констант являются целыми числами, и у вас есть куча 1/<constant> целочисленных делений в вашем уравнении.

Действительно, если я заменю / на // в моей версии (для Python 3)Я могу воспроизвести вашу проблему.

Просто добавьте from __future__ import division вверху вашего скрипта, чтобы решить вашу проблему.

Также добавьте from __future__ import print_function вв верхней части замените print <something> на print(<something>), и ваш сценарий полностью совместим с Python 3 и 2).

...