#numpy #scipy #ode #runge-kutta
Вопрос:
Я хотел бы решить уравнения ОДУ первого порядка движения с помощью функции scipy solve_ivp. Я вижу, что делаю что-то не так, потому что это должен быть эллипс, но я рисую только четыре точки. Способны ли вы распознать ошибку?
import math import matplotlib.pyplot as plt import numpy as np import scipy.integrate gim = 4*(math.pi**2) x0 = 1 #x-position of the center or h y0 = 0 #y-position of the center or k vx0 = 0 #vx position vy0 = 1.1* 2* math.pi #vy position initial = [x0, y0, vx0, vy0] #initial state of the system time = np.arange(0, 1000, 0.01) #period def motion(t, Z): dx = Z[2] # vx dy = Z[3] # vy dvx = -gim/(x**2 y**2)**(3/2) * x * Z[2] dvy = -gim/(x**2 y**2)**(3/2) * y * Z[3] return [dx, dy, dvx, dvy] sol = scipy.integrate.solve_ivp(motion, t_span=time, y0= initial, method='RK45') plt.plot(sol.y[0],sol.y[1],"x", label="Scipy RK45 solution") plt.show() 
Комментарии:
1. Я не думаю, что вы опубликовали точный код, который вы используете:
xиyне определены вmotion.solve_ivpаргументt_spanожидает конечные точки временного интервала (т. е. два числа). Наконец, что это за моделирование? Наличие зависимости dvx от vx (и аналогично для dvy) приводит к затуханию, которое, я не думаю , приведет к эллипсу; однако уравнения нелинейны, поэтому трудно сказать.
Ответ №1:
Код не должен запускаться из новой рабочей области. Переменные x,y , которые используются в формуле гравитации, нигде не объявляются. Поэтому вставьте строку x,y = Z[:2] или что-то подобное.
Формула гравитации обычно не содержит компонент скорости. Удалить Z[2] , Z[3] .
Проверьте еще раз, какой промежуток времени и время оценки вы ожидаете. Временной интервал принимает первые два значения из массива. Поэтому измените значение t_span=time[[0,-1]] на создание пары первого и последнего значения времени.
Второй участок страдает от недостаточного количества оценочных точек, используемые отрезки слишком велики. С вашим time массивом это не должно быть проблемой.