Метод Эйлера в С++; Значения становятся слишком большими слишком быстро

Я пытаюсь решить уравнение движения для частицы с массой 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 итераций, и я не могу понять, почему. Когда я рисую позицию как функцию времени, я получаю график ниже, что явно неверно.

Ваше умножение на t, вероятно, должно быть на delta вместо этого

Mat 28.10.2022 15:25

Вы должны изменить свой код, чтобы использовать непосредственно тестовые данные, или предоставить тестовые данные, чтобы можно было дублировать ваши наблюдения. // Иногда помогает "интегрировать" время как дополнительный компонент с производной константой 1. Если есть явные несоответствия, шаг метода имеет некоторую ошибку реализации.

Lutz Lehmann 28.10.2022 15:32

Вы уверены в условии (x_new != x_new) || (v_new != v_new) ??

Yves Daoust 02.11.2022 23:42

Вы забыли сообщить нам входные параметры.

Yves Daoust 02.11.2022 23:46
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
2
4
73
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Ваша реализация уравнения неверна. Вы используете 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.

Yves Daoust 02.11.2022 23:50

Другие вопросы по теме