Интерпретация вывода dsolve для использования с ODE45 [MATLAB]

Question

Я пытался решить дифференциальное уравнение второго порядка, используя это 1, , однако я не мог понять это правильно и не нашел ничего полезного в Интернете, но я считаю, что я достиг прогресса.

Я использовал dsolve;

syms x(t) v(t) fi(t) 

[x(t), v(t)] = dsolve(diff(x) == v, diff(v) == fi/m, x(0) == [-L, -L], v(0) == [5, 10]) 

который дает мне;

x(t) = 



int(fi(x)/5, x, 0, t, 'IgnoreSpecialCases', true, 'IgnoreAnalyticConstraints', true) - 5 





v(t) = 

C2 + t*(int(fi(x)/5, x, 0, t, 'IgnoreSpecialCases', true, 'IgnoreAnalyticConstraints', true) - 5) + int(-(x*fi(x))/5, x, 0, t, 'IgnoreSpecialCases', true, 'IgnoreAnalyticConstraints', true) 

Теперь мне нужна помощь, связанная с результатом, и мне интересно, могу ли я использовать этот результат, чтобы получить что-то из ode45? Также я хочу построить решение в качестве эталонной траектории для моделирования 500 частиц, движущихся через силовое поле.

---

использованием ode45:

function dxdt = solution(t,y0) 
frprintf('Second stop') 
..... 
dxdt = [x, v] 
end 

вызова из основного файла:

t:dt:t_f 
y0 = [x0,v0] 
fprintf('first stop') 
[x, v] = ode45(@solution, y0, t) 

Я настроил его так, что если код будет работать гладко было бы напечатать «первую остановку, вторая остановка, третья остановка "и" четвертая остановка ", она только распечатывает первую остановку и вот где я получаю ошибку.

Answer 1

Предполагая, что вы где-то определено

function F = fi(t) 
    F = ... 
end 

вы определяете функцию ОДУ как (с помощью m в качестве глобальной переменной)

function doty = odefunc(t,y) 
    doty = [ y(2); fi(t)/m ] 
end 

и вызвать

t = t0:dt:tf 
y0 = [ x0, v0 ] 
t,y = ode45(odefunc, t, y0) 

plot(t, y(:,1)) 

---

Как правило, V ector y будет содержать точку в фазовом пространстве (положение и скорость/импульс) для всех частиц или объектов в системе. Затем в odefunc вы вычисляете силы для этой конкретной фазовой пространственной точки в момент времени t и составляете из нее вектор производной в точку фазового пространства.

В 3D-моделировании вы, например, могли бы организовать, что y(6*(k-1)+1:6*(k-1)+3) - это положение частицы k и y(6*(k-1)+4:6*(k-1)+6) вектор скорости. Или вы можете разделить положение и скорость с помощью y(3*(k-1)+1:3*(k-1)+3) для положения и y(3*(N+k-1)+1:3*(N+k-1)+3) для скорости для k-ой из N частиц.