Я пытаюсь решить уравнение движения для частицы с массой m, прикрепленной к пружине с коэффициентом пружины k. Однако оба установлены на 1. Алгоритм выглядит так:
Мое (попытка) решение, написанное на С++, выглядит так:
#include <iostream>
#include <iomanip>
#include <math.h>
#include <stdlib.h>
#include <fstream>
// Initialise file to write series of values in
std::ofstream output("Eulermethod.txt");
// Define Euler algorithm
void euler(double x_0, double v_0, double delta, double t_max) {
double x_prev = x_0;
double v_prev = v_0;
double x_new, v_new;
for (double t = 0; t < t_max; t = t + delta) {
x_new = x_prev + t * v_prev;
v_new = v_prev - t * x_prev;
// Writes time, position and velocity into a csv file
output << std::fixed << std::setprecision(3) << t << "," << x_prev << "," << v_prev << std::endl;
x_prev = x_new;
v_prev = v_new;
// Breaks loop if values get to big
if ((x_new != x_new) || (v_new != v_new) || (std::isinf(x_new) == true) || (std::isinf(v_new) == true)) {
break;
}
}
}
int main() {
// Initialize with user input
double x_0, v_0, t_max, delta;
std::cout << "Initial position x0?: ";
std::cin >> x_0;
std::cout << "Intial velocity v0?: ";
std::cin >> v_0;
std::cout << "Up to what time t_max?: ";
std::cin >> t_max;
std::cout << "Step size delta?: ";
std::cin >> delta;
// Runs the function
euler(x_0, v_0, delta, t_max);
}
Я знаю, что решение будет расти бесконечно, но для меньших значений t оно должно напоминать аналитическое решение, но расти медленно. Значения, которые я получаю, выходят из пропорций после ок. 10 итераций, и я не могу понять, почему. Когда я рисую позицию как функцию времени, я получаю график ниже, что явно неверно.
Вы должны изменить свой код, чтобы использовать непосредственно тестовые данные, или предоставить тестовые данные, чтобы можно было дублировать ваши наблюдения. // Иногда помогает "интегрировать" время как дополнительный компонент с производной константой 1. Если есть явные несоответствия, шаг метода имеет некоторую ошибку реализации.
Вы уверены в условии (x_new != x_new) || (v_new != v_new) ??
Вы забыли сообщить нам входные параметры.
Ваша реализация уравнения неверна. Вы используете t
вместо dt
. Правильный вариант:
x_new = x_prev + delta * v_prev;
v_new = v_prev - delta * x_prev;
И примечание, если вы планируете развивать свой код дальше: общий подход к реализации решателя ODE заключается в том, чтобы иметь метод с сигнатурой, подобной
Output = solveOde(System, y0, t);
Где System
— метод, описывающий ОДУ dy/dx = f(x,t)
, например.
std::vector<double> yourSystem(std::vector<double> y, double /*t unused*/)
{
return {y[1], -y[0]};
}
y0
— начальные условия, а t
— вектор времени (дельта вычисляется внутри). Взгляните на boost odeint или более компактную и прозрачную документацию по python.
В реализованном виде этот код решает tx" - x' + t³x = 0 и находит тривиальное решение x = 0.
Ваше умножение на
t
, вероятно, должно быть наdelta
вместо этого