from numpy import*
import matplotlib.pyplot as plt
from scipy.integrate import *
def odein():
#dy1/dt=y2
#dy2/dt=y1**2+1:
def f(y,t):
return y**2+1
t =arange(0,1,0.01)
y0 =0.0
y=odeint(f, y0,t)
y = array(y).flatten()
return y,t
def oden():
f = lambda t, y: y**2+1
ODE=ode(f)
ODE.set_integrator('dopri5')
ODE.set_initial_value(0, 0)
t=arange(0,1,0.01)
z=[]
t=arange(0,1,0.01)
for i in arange(0,1,0.01):
ODE.integrate(i)
q=ODE.y
z.append(q[0])
return z,t
def rungeKutta(f, to, yo, tEnd, tau):
def increment(f, t, y, tau):
if z==1:
k0 =tau* f(t,y)
k1 =tau* f(t+tau/2.,y+k0/2.)
k2 =tau* f(t+tau/2.,y+k1/2.)
k3 =tau* f(t+tau, y + k2)
return (k0 + 2.*k1 + 2.*k2 + k3) / 6.
elif z==0:
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6
t = []
y= []
t.append(to)
y.append(yo)
while to < tEnd:
tau = min(tau, tEnd - to)
yo = yo + increment(f, to, yo, tau)
to = to + tau
t.append(to)
y.append(yo)
return array(t), array(y)
def f(t, y):
f = zeros([1])
f[0] = y[0]**2+1
return f
to = 0.
tEnd = 1
yo = array([0.])
tau = 0.01
z=1
t, yn = rungeKutta(f, to, yo, tEnd, tau)
y1n=[i[0] for i in yn]
plt.figure()
plt.title("Абсолютная погрешность численного решения(т.р.- u(t)=tan(t)) ДУ\ndu/dt=u**2+1 c u(0)=0 при t>0")
plt.plot(t,abs(array(y1n)-array(tan(t))),label='Метод Рунге—Кутта \nчетвертого порядка - расчёт по алгоритму')
plt.xlabel('Время')
plt.ylabel('Абсолютная погрешность.')
plt.legend(loc='best')
plt.grid(True)
z=0
t, ym = rungeKutta(f, to, yo, tEnd, tau)
y1m=[i[0] for i in ym]
plt.figure()
plt.title("Абсолютная погрешность численного решения(т.р.- u(t)=tan(t)) ДУ\ndu/dt=u**2+1 c u(0)=0 при t>0")
plt.plot(t,abs(array(y1m)-array(tan(t))),label='Метод Рунге—Кутта— Фельберга \nпятого порядка - расчёт по алгоритму')
plt.xlabel('Время')
plt.ylabel('Абсолютная погрешность.')
plt.legend(loc='best')
plt.grid(True)
plt.figure()
plt.title("Абсолютная погрешность численного решения (т.р.- u(t)=tan(t)) ДУ\ndu/dt=u**2+1 c u(0)=0 при t>0")
y,t=odein()
plt.plot(t,abs(array(tan(t))-array(y)),label='Функция odein')
plt.xlabel('Время')
plt.ylabel('Абсолютная погрешность.')
plt.legend(loc='best')
plt.grid(True)
plt.figure()
plt.title("Абсолютная погрешность численного решения (т.р.- u(t)=tan(t)) ДУ\ndu/dt=u**2+1 c u(0)=0 при t>0")
z,t=oden()
plt.plot(t,abs(tan(t)-z),label='Функция ode метод Рунге—Кутта— Фельберга \nпятого порядка')
plt.xlabel('Время')
plt.ylabel('Абсолютная погрешность.')
plt.legend(loc='best')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import time
start = time.time()
from scipy.integrate import ode
ts = [ ]
ys = [ ]
FlightTime, Distance, Height =0,0,0
y4old=0
def fout(t, y):# обработчик шага
global FlightTime, Distance, Height,y4old
ts.append(t)
ys.append(list(y.copy()))
y1, y2, y3, y4 = y
if y4*y4old<=0: # достигнута точка максимума
Height=y3
if y4<0 and y3<=0.0: # тело достигло поверхности
FlightTime=t
Distance=y1
return -1
y4old=y4
# функция правых частей системы ОДУ
def f(t, y, k): # имеется дополнительный аргумент k
g=9.81
y1, y2, y3, y4 = y
return [y2,-k*y2*np.sqrt(y2**2+y4**2), y4,-k*y4*np.sqrt(y2**2+y4**2)-g]
tmax=1.41 # максимально допустимый момент времени
alph=np.pi/4 # угол бросания тела
v0=10.0 # начальная скорость
K=[0.1,0.2,0.3,0.5] # анализируемые коэффициенты сопротивления
y0,t0=[0, v0*np.cos(alph), 0, v0*np.sin(alph)], 0 # начальные условия
ODE=ode(f)
ODE.set_integrator('dopri5', max_step=0.01)
ODE.set_solout(fout)
fig, ax = plt.subplots()
fig.set_facecolor('white')
for k in K: # перебор значений коэффициента сопротивления
ts, ys = [ ],[ ]
ODE.set_initial_value(y0, t0) # задание начальных значений
ODE.set_f_params(k) # передача дополнительного аргумента k # в функцию f(t,y,k) правых частей системы ОДУ
ODE.integrate(tmax) # решение ОДУ
print('Flight time = %.4f Distance = %.4f Height =%.4f '% (FlightTime,Distance,Height))
Y=np.array(ys)
plt.plot(Y[:,0],Y[:,2],linewidth=3,label='k=%.1f'% k)
stop = time.time()
plt.title("Результаты численного решения системы четырёх ОДУ \n с использованием функци ode с атрибутом dopri5 ")
print ("Время на модельную задачу: %f"%(stop-start))
plt.grid(True)
plt.xlim(0,8)
plt.ylim(-0.1,2)
plt.legend(loc='best')
plt.show()
def increment(f, t, y, tau
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6
from numpy import*
import matplotlib.pyplot as plt
import time
start = time.time()
def rungeKutta(f, to, yo, tEnd, tau):
def increment(f, t, y, tau):# поиск приближённого решения методом Рунге—Кутта—Фельберга.
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6
t = []#подготовка пустого списка t
y= []#подготовка пустого списка y
t.append(to)#внесение в список t начального значения to
y.append(yo)#внесение в список y начального значения yo
while to < tEnd:#внесение результатов расчёта в массивы t,y
tau = min(tau, tEnd - to)#определение минимального шага tau
yo = yo + increment(f, to, yo, tau) # расчёт значения в точке t0,y0 для задачи Коши
to = to + tau # приращение времени
t.append(to) # заполнение массива t
y.append(yo) # заполнение массива y
return array(t), array(y)
def f(t, y): # функция правых частей системы ОДУ
f = zeros([4])
f[0]=y[1]
f[1]=-k*y[1]*sqrt(y[1]**2+y[3]**2)
f[2]=y[3]
f[3]=-k*y[3]*sqrt(y[1]**2+y[3]**2) -g
if y[3]<0 and y[2]<=0.0: # тело достигло поверхности
return -1
return f
to = 0# начальный момент отсчёта времени
tEnd = 1.41# конечный момент отсчёта времени
alph=pi/4# угол бросания тела
v0=10.0 # начальная скорость
K=[0.1,0.2,0.3,0.5]# анализируемые коэффициенты сопротивления среды
g=9.81
yo = array([0.,v0*cos(alph),0.,v0*sin(alph)]) # начальные условия
tau =0.01# шаг
for i in K: # перебор значений коэффициента сопротивления среды
k=i
t, y = rungeKutta(f, to, yo, tEnd, tau)
y1=array([i[0] for i in y]) # извлечение переменных из массива y
y3=array([i[2] for i in y])
# визуализация сопротивление среды "к" и расчёт и визуализация s,h,t
plt.plot(y1,y3,linewidth=2,label='k=%.1f h=%.3f s=%.2f t=%s' % (k,max(y3),max(y1),round(t[list(y1).index(max(y1))],3)))
stop = time.time()
plt.title("Результаты численного решения системы четырёх ОДУ \n с использованием адаптированного для Python\n метода Рунге—Кутта—Фельберга ")
print ("Время на модельную задачу: %f"%(stop-start))
plt.xlabel('Высота h')
plt.ylabel('Расстояние s')
plt.legend(loc='best')
plt.xlim(0,8)
plt.ylim(-0.1,2)
plt.grid(True)
plt.show()
# Метод пристрелки
from numpy import*
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm,os
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
from scipy.integrate import odeint
from scipy import linalg
import time
start = time.time()
c1 = 1.0
c2 = 0.8
c3 = 0.5
a =0.0
b = 1.0
nn =100
initial_state_0 =array( [a, c1, 0.0, 0.0])
initial_state_I =array( [a, 0.0, 1.0, 0.0])
initial_state_II =array( [a, 0.0, 0.0, 1.0])
to = a
tEnd =b
N = int(nn)
tau=(b-a)/N
def rungeKutta(f, to, yo, tEnd, tau):
def increment(f, t, y, tau):
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6
t = []
y= []
t.append(to)
y.append(yo)
while to < tEnd:
tau = min(tau, tEnd - to)
yo = yo + increment(f, to, yo, tau)
to = to + tau
t.append(to)
y.append(yo)
return array(t), array(y)
def f(t, y):
global theta
f = zeros([4])
f[0] = 1
f[1] = -y [1]-y[2] +theta* sin(y[0])
f[2] = -y[2]+y[3]
f[3] = -y[2]
return f
# Решение НЕОДНОРОДНОЙ системы -- theta = 1
theta = 1.0
yo =initial_state_0
t, y = rungeKutta(f, to, yo, tEnd, tau)
y2=[i[2] for i in y]
y3=[i[3] for i in y]
# Извлечение требуемых для решения задачи значений
# Y20 = Y2(b), Y30 = Y3(b)
Y20 = y2[N-1]
Y30 = y3[N-1]
# Решение ОДНОРОДНОЙ системы -- theta = 0, задача I
theta = 0.0
yo= initial_state_I
t, y = rungeKutta(f, to, yo, tEnd, tau)
y2=[i[2] for i in y]
y3=[i[3] for i in y]
# Извлечение требуемых для решения задачи значений
# Y21= Y2(b), Y31 = Y3(b)
Y21= y2[N-1]
Y31 = y3[N-1]
# Решение ОДНОРОДНОЙ системы -- theta = 0, задача II
theta = 0.0
yo =initial_state_II
t, y = rungeKutta(f, to, yo, tEnd, tau)
y2=[i[2] for i in y]
y3=[i[3] for i in y]
# Извлечение требуемых для решения задачи значений
# Y211= Y2(b), Y311 = Y3(b)
Y211= y2[N-1]
Y311 = y3[N-1]
# Формирование системы линейных
# АЛГЕБРАИЧЕСКИХ уравнений для определния p2, p3
b1 = c2 - Y20
b2 = c3 - Y30
A = array([[Y21, Y211], [Y31, Y311]])
bb = array([[b1], [b2]])
# Решение системы
p2, p3 = linalg.solve(A, bb)
# Окончательное решение краевой
# НЕОДНОРОДНОЙ задачи, theta = 1
theta = 1.0
yo = array([a, c1, p2, p3])
t, y = rungeKutta(f, to, yo, tEnd, tau)
y0=[i[0] for i in y]
y1=[i[1] for i in y]
y2=[i[2] for i in y]
y3=[i[3] for i in y]
# Проверка
print('y0[0]=', y0[0])
print('y1[0]=', y1[0])
print('y2[0]=', y2[0])
print('y3[0]=', y3[0])
print('y0[N-1]=', y0[N-1])
print('y1[N-1]=', y1[N-1])
print('y2[N-1]=', y2[N-1])
print('y3[N-1]=', y3[N-1])
j = N
xx = y0[:j]
yy1 = y1[:j]
yy2 = y2[:j]
yy3 = y3[:j]
stop = time.time()
print ("Время на модельную задачу: %f"%(stop-start))
plt.subplot(2, 1, 1)
plt.plot([a], [c1], 'ro')
plt.plot([b], [c2], 'go')
plt.plot([b], [c3], 'bo')
plt.plot(xx, yy1, color='r') #Построение графика
plt.plot(xx, yy2, color='g') #Построение графика
plt.plot(xx, yy3, color='b') #Построение графика
plt.xlabel(r'$x$') #Метка по оси x в формате TeX
plt.ylabel(r'$y_k(x)$') #Метка по оси y в формате TeX
plt.title(r'Метод пристрелки ', color='blue')
plt.grid(True) #Сетка
patch_y1 = mpatches.Patch(color='red',
label='$y_1$')
patch_y2 = mpatches.Patch(color='green',
label='$y_2$')
patch_y3 = mpatches.Patch(color='blue',
label='$y_3$')
plt.legend(handles=[patch_y1, patch_y2, patch_y3])
ymin, ymax = plt.ylim()
xmin, xmax = plt.xlim()
plt.subplot(2, 1, 2)
font = {'family': 'serif',
'color': 'blue',
'weight': 'normal',
'size': 12,
}
plt.text(0.2, 0.8, r'$\frac{dy_1}{dx}= - y_1 - y_2 + \sin(x),$',
fontdict=font)
plt.text(0.2, 0.6,r'$\frac{dy_2}{dx}= - y_1 + y_3,$',
fontdict=font)
plt.text(0.2, 0.4, r'$\frac{dy_3}{dx}= - y_2 - y_2,$',
fontdict=font)
plt.text(0.2, 0.2, r'$y_1(a)=c_1, '
r'\quad y_2(b)=c_2, \quad y_3(b)=c_3.$',
fontdict=font)
plt.subplots_adjust(left=0.15)
plt.show()
К сожалению, не доступен сервер mySQL