Рунге-Кутта ведет себя странно, когда начальное значение отрицательное

У меня есть следующий алгоритм RK4 на C++ для решения дифференциальных уравнений первого порядка:

using ODE_Function = std::function<double/*dy/dx*/(double/*x*/, double/*y*/)>;

template<int length> //"length" is the number of data points
double* rk4(ODE_Function fxn, double y0, double x0, double h) {

    double* y = new double[length];
    
    y[0] = y0; //initialize output array
    double x = x0; //n is our discretized x variable

    //main loop
    for (int n = 1; n < length; n++) {//y0 is already known, so we start from the second data point
        
        double k1 = fxn(x, y[n-1]);
        double k2 = fxn(x + h/2, y[n-1] + h*k1/2);
        double k3 = fxn(x + h/2, y[n-1] + h*k2/2);
        double k4 = fxn(x + h, y[n-1] + h*k3);

        y[n] = y[n-1] + h/6*(k1 + 2*k2 + 2*k3 + k4);
        x += h;
    }

    return y;
}

Однако это ведет себя странно, когда x0 отрицательно...

В частности, следующие случаи теоретически должны выводить результаты, которые следуют тенденции (-x^2)-2(x+1), просто просматриваемой из разных «окна». Но это делает только финальный. Первый приводит к отрицательной экспоненте, а второй по какой-то причине приводит к положительной экспоненте.

double fxn(double x, double y) {
    return x*x + y;
}
int main() {
    const int length = 100;
    double y0 = -2;
    double x0 = -1; //LOOK HERE
    double h = 0.1;

    double* result = rk4<length>(fxn, y0, x0, h);
}
double fxn(double x, double y) {
    return x*x + y;
}
int main() {
    const int length = 100;
    double y0 = -2;
    double x0 = -5; //LOOK HERE
    double h = 0.1;

    double* result = rk4<length>(fxn, y0, x0, h);
}
double fxn(double x, double y) {
    return x*x + y;
}
int main() {
    const int length = 100;
    double y0 = -2;
    double x0 = 0; //LOOK HERE
    double h = 0.1;

    double* result = rk4<length>(fxn, y0, x0, h);
}

какой у Вас вопрос?

463035818_is_not_an_ai 25.06.2024 10:20

почему-то... -- Что такое отладчик?. Вам нужно выполнить расчет вручную, а затем сравнить результаты при отладке вашей программы, чтобы увидеть, где (и когда) ваша реализация отличается от того, что вы рассчитали вручную.

PaulMcKenzie 25.06.2024 11:22
Стоит ли изучать 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 называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
0
2
79
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

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

В вашей реализации Рунге-Кутты нет ничего плохого — просто ваше ожидание решения.

Немного переформулируйте уравнение и используйте интегрирующий коэффициент. Вы найдете точное решение

y = (y0+x02+2x0+2)exp(x-x0)-x2-2x-2

Попробуйте это с помощью функции main(). (Вам нужно будет включить заголовки iostream, function и cmath.)

int main() {
    const int length = 100;
    double y0 = -2;
    double x0 = -5; //LOOK HERE
    double h = 0.1;

    double* result = rk4<length>(fxn, y0, x0, h);

    for ( int i = 0; i < length; i++ ) 
    {
       double x = x0 + i * h;
       double exact = ( y0 + x0 * x0 + 2 * x0 + 2 ) * exp( x - x0 ) - x * x - 2 * x - 2;
       cout << i << '\t' << x << '\t' << result[i] << '\t' << exact << '\n';
    }
}

Большое спасибо за ваш ответ! Однако когда я вычисляю «ожидаемую» функцию, перед экспонентой стоит константа, которая определяется начальным условием y(0). Когда y(0)=-2, эта константа должна быть равна 0, оставляя только квадратичный член решения, но мой код Рунге-Кутты не делает это правильно. Для справки обратите внимание, что в: wolframalpha.com/… c1 = 0 для y(0)= -2

Aryan MP 25.06.2024 18:59

Вы можете проверить простым дифференцированием, что это решает dy/dx=x^2+y и проходит через точку (x0,y0). Ваша константа c1 равна y0+x0^2+2x0+2. Это явно зависит от x0. В противном случае вы бы не достигли граничного условия.

lastchance 25.06.2024 19:45

Если вы установите x0=-5, то вы применяете начальное граничное условие y(-5)=-2, а НЕ y(0)=-2. Так что ваш сценарий WolframAlpha не имеет значения.

lastchance 25.06.2024 20:48

Оооооооооооооооооооооооооооооо ты так прав

Aryan MP 25.06.2024 20:52

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