Появилась потребность попробовать в деле scipy и numpy, а заодно и matplotlib в задачах вычислительной физики. Взял за пример задачу решить уравнение Шрёдингера для одномерной квантовой частицы в потенциальной яме конечной высоты . Забегая вперёд, скажу, эти python-библиотеки справились с решением весьма быстро и качественно. Получившийся код, при надлежащей модификации, вполне можно использовать для решения целого класса дифференциальных уравнений второго порядка.
Итак, после очевидных трансформаций, имеем одномерное уранение Шрёдингера:
Для его решения будем использовать интегрирующую функцию odeint из модуля scipy.integrate . Однако для этого нам понадобится преобразовать данное уравнение в систему дифференциальных уравнений первого порядка. Опустим нехитрые математические выкладки (которые можно провести самостоятельно) и получим в итоге систему уравнений первого порядка:
при
С получением данной системы уравнений математическая часть задачи закончена, можно начинать кодировать.
Python — численное решение дифференциального уравнения 1го порядка и вывод графика
Для начала, импортируем всё что нужно
from scipy import integrate import numpy as np # на будущее import math # . import matplotlib.pyplot as plt
Затем, закодируем переменные задачи
# энергия частицы e=1.3 # константа частицы h=0.5 # фактор потенциальной ямы A=1.5 # константы интегрирования a = 0.0 b = 10.0 step = 0.005 # начальные условия alpha=0 beta=1 # неиспользуемые дополнительные переменные omega=1 phi=0
Пропишем функцию потенциала
def u(x): «»» Функция потенциала, прямоугольная потенциальная яма «»» if x>5.0: return A*1e1 else: return A
Функцию, задающую систему уравнений
def se_solve(Y,x): «»» Функция, задающая систему дифференциальных уравнений «»» return [Y[1],(e-u(x))/h*Y[0]]
def main(): # массивы для решения задачи a_x=np.arange(a,b,step) y=np.arange(a,b,step) # интегрируем систему уравнений asol=integrate.odeint(se_solve,[alpha,beta],a_x) # на графике будем изображать квадрат модуля волновой функции # несущий физический смысл плотности вероятности for i in range(0,len(asol)): y[i]=asol[i][0]**2 plt.plot(y) plt.show() if __name__ == ‘__main__’: main()
Запускаем и получаем ожидаемый результат: график волновой функции квантовой частицы в одномерной потенциальной яме конечной высоты с туннельным эффектом.
Исходник доступен на гитхабе, всем добра.
Источник: flaurion.livejournal.com
Численное решение обыкновенных дифференциальных уравнений (ОДУ) в Python
Модуль scipy.integrate имеет две функции ode() и odeint(), которые предназначены для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (т.е. задача Коши).
Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.
Python для самых маленьких. Линейные уравнения. Решение задач
from scipy.integrate import odeint
Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат
odeint(func, y0, t[,args=(), . ])
Решение одного ОДУ
Допустим надо решить диф. уравнение 1-го порядка
import numpy as np from scipy. integrate import odeint import matplotlib.pyplot as plt # create function def dydt(y, t): return -y*t t = np.linspace( -2, 2, 51) # vector of time y0 = 1 # start value y = odeint (dydt, y0, t) # solve eq. y = np.array(y).flatten() plt.plot( t, y,’-sr’, linewidth=3) # graphic plt.show() # display
Получилось что-то такое:
Решение системы ОДУ
Пусть теперь мы хотим решить (автономную) систему диф. уравнений 1-го порядка
import numpy as np from scipy. integrate import odeint import matplotlib.pyplot as plt # create function def f(y, t): y1, y2 = y return [y2, — y2 — y1] t = np.linspace( 0, 10, 41) # vector of time y0 = [0, 1] # start value w = odeint(f, y0, t) # solve eq. y1 = w[:,0] y2 = w[:,1] fig = plt.figure(facecolor=’white’) plt.plot(t, y1, ‘-o’, t, y2, ‘-o’, linewidth=2) plt.ylabel(«z») plt.xlabel(«t») plt.grid(True) plt.show() # display
Выходной массив w состоит из двух столбцов — y1(t) и y2(t).
Также без труда можно построить фазовые траектории:
fig2 = plt.figure(facecolor=’white’) plt.plot(y1, y2, linewidth=2) plt.grid(True) plt.show()
Платы ARDUINO
Now 05.07.23 6:04:08, Your IP: 46.158.74.212; spyphy.zl3p.com/python/34_scipy_ode
Источник: spyphy.zl3p.com
Русские Блоги
Библиотека sympy для поиска дифференциальных уравнений в Python
Вопрос 1:
f ′ ′ ( x ) − 2 f ′ ( x ) + f ( x ) = s i n ( x ) f»(x)-2f'(x) + f(x) = sin(x) f ′ ′ ( x ) − 2 f ′ ( x ) + f ( x ) = s i n ( x )
Порядок действий следующий
from sympy import * f= symbols(‘f’, cls=Function) eq = Eq(f(x).diff(x, x) — 2*f(x).diff(x) + f(x), sin(x)) print(dsolve(eq, f(x)))
Eq(f(x), (C1 + C2*x)*exp(x) + cos(x)/2)
~
~
~
~
Вложение:Назначьте два вопроса на экзамене
1. Используйте библиотеку Sympy Python для решения дифференциальных уравнений. y = f ( x ) y=f(x) y = f ( x ) И попробуйте использовать matplotlib для рисования изображений функций
f ′ ( x ) + f ( x ) + f 2 ( x ) = 0 , f ( 0 ) = 1 f'(x)+f(x)+f^2(x)=0,qquad f(0)=1 f ′ ( x ) + f ( x ) + f 2 ( x ) = 0 , f ( 0 ) = 1
Порядок действий следующий
from sympy import * f = symbols(‘f’, cls=Function) x = symbols(‘x’) eq = Eq(f(x).diff(x,1)+f(x)+f(x)**2, 0) print(dsolve(eq, f(x))) C1 = symbols(‘C1’) eqr = -C1/(C1 — exp(x)) eqr1 = eqr.subs(x, 0) print(solveset(eqr1 — 1, C1)) eqr2 = eqr.subs(C1, 1/2) # Рисовать import matplotlib.pyplot as plt import numpy as np x_1 = np.arange(-5, 5, 0.1) y_1 = [-0.5/(0.5 — exp(x)) for x in x_1] plt.plot(x_1, y_1) plt.axis([-6,6,-10,10]) plt.grid() plt.show()
Eq(f(x), -C1/(C1 — exp(x))) FiniteSet(1/2)
~
2. Используйте библиотеку Sympy Python для решения дифференциальных уравнений. y = y ( x ) y=y(x) y = y ( x ) И попробуйте использовать matplotlib для рисования изображений функций
y ′ ( x ) = y ( x ) , y ( 0 ) = 1 y'(x)=y(x),qquad y(0)=1 y ′ ( x ) = y ( x ) , y ( 0 ) = 1
Порядок действий следующий
from sympy import * y = symbols(‘y’, cls=Function) x = symbols(‘x’) eq = Eq(y(x).diff(x,1), y(x)) print(dsolve(eq, y(x))) C1 = symbols(‘C1′) eqr = C1*exp(x) eqr1 = eqr.subs(x, 0) print(solveset(eqr1 — 1, C1)) eqr2 = eqr.subs(C1, 1) # Рисовать import matplotlib.pyplot as plt import numpy as np x_1 = np.arange(-5, 5, 0.01) y_1 = [exp(x) for x in x_1] plt.plot(x_1, y_1, color=’orange’) plt.grid() plt.show()
Eq(y(x), C1*exp(x)) FiniteSet(1)
Источник: russianblogs.com