Я пытаюсь выполнить уравнения Лоренца в питоне (я следую упражнению 8.3 от Марка Ньюмана - Вычислительная физика (2012, Независимая издательская платформа CreateSpace)) Я уже получил графику, и все выглядит "правильно". Это, вероятно, математическая проблема, а не проблема программирования, но я пишу здесь, чтобы убедиться. Во-первых, это мой код:
from numpy import arange,array
import pylab
def f(v,t):
s=10
r=28
b=8/3
x= v[0]
y= v[1]
z= v[2]
fx= s*(y - x)
fy= r*x - y - x*z
fz= x*y - b*z
return array([fx,fy,fz],float)
def d(N):
a=0.0
b=50.0
h=(b-a)/N
r=array([0.0,1.0,0.0],float)
tpoints=arange(a,b,h)
xpoints= []
ypoints= []
zpoints= []
for t in tpoints:
xpoints.append(r[0])
ypoints.append(r[1])
zpoints.append(r[2])
k1 = h*f(r,t)
k2 = h*f(r+0.5*k1,t+0.5*h)
k3 = h*f(r+0.5*k2,t+0.5*h)
k4 = h*f(r+k3,t+h)
r += (k1+2*k2+2*k3+k4)*(1/6)
return tpoints,xpoints,ypoints,zpoints
for i in range (1,6):
N=10**i
pylab.plot(d(N)[0],d(N)[1], label=N)
pylab.xlabel("t")
pylab.ylabel("x(t)")
pylab.title("Gráficos x em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[0],d(N)[2], label=N)
pylab.xlabel("t")
pylab.ylabel("y(t)")
pylab.title("Gráficos y em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[0],d(N)[3], label=N)
pylab.xlabel("t")
pylab.ylabel("z(t)")
pylab.title("Gráficos z em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[1],d(N)[3], label=N)
pylab.xlabel("x")
pylab.ylabel("z(x)")
pylab.title("Gráficos z em função de x")
pylab.legend()
pylab.show()
Это дает мне графики для решения проблемы, и я думаю, что это правильно.
Когда я перехожу от i=1
к i=3
, я получаю эту ошибку:
Я думаю, что это связано с математической проблемой, но когда я ищу ошибку, я получаю что-то с массивами. Так что я проверяю это.
Для i
равного или более 4 это не доставило мне проблем.
1- Я думаю, что моя проблема в том, что я понимаю, что означают эти предупреждения. Я даже подумал, что это ошибки. 2-что вы имеете в виду, когда ставите "точку останова прямо перед каждым предупреждением" 3-"каковы все значения в этих уравнениях?" не понял этого. 4- Если вы имеете в виду, что я получаю то, что хочу, с N => 10000 да, с N <= 1000, нет.
Вы программируете и не знаете о точках останова? Обучение использованию отладчика похоже на важную информацию первого дня.
@RandomDavis: если при использовании jupyter нет точек останова, потребуется подходящая IDE, такая как pyzo или spyder и т. д., чтобы иметь обычные средства отладки. В jupyter можно было вставлять операторы печати внутри циклов, чтобы увидеть список до последних вычисленных значений.
@LutzLehmann import pdb; pdb.set_trace()
работает в любом переводчике
@PaulH: Спасибо, я это запомню.
Я использую jupyter, но это все равно математическая «проблема», спасибо.
Система Лоренца с RK4 требует размера шага около 0.05
или меньше, то есть N=10**4
или больше для вашей конструкции. Смотрите почти дубликат аттрактора Лоренца с питоном Рунге-Кутта.
Для больших размеров шага, то есть случаев, в которых вы получаете ошибку, метод будет возвращать хаотические результаты, которые в основном не связаны с точным решением системы и любыми связанными с ним границами. Таким образом, расхождение с переполнением с плавающей запятой возможно из-за квадратичных суперлинейных членов.
Вы понимаете, что означают эти предупреждения? Кроме того, если вы поставите точку останова прямо перед каждым предупреждением, каковы все значения в этих уравнениях? Они то, что вы ожидаете?