математика — алгоритм Рунге-Кутты переполнение стека

Question

математика — алгоритм Рунге-Кутты переполнение стека

Ниже приведен мой алгоритм Рунге-Кутты 4-го порядка для решения ОДУ первого порядка. Я проверяю это по найденному примеру википедии Вот решать:

\frac{dx}{dt} = tan(x) + 1

К сожалению это немного. Я долго играл, но не могу найти ошибку. Ответ должен быть t = 1.1 и x = 1.33786352224364362. Приведенный ниже код дает t = 1,1 и x = 1,42223.

/*

This code is a 1D classical Runge-Kutta method. Compare to the Wikipedia page.

*/

#include <math.h>
#include <iostream>
#include <iomanip>

double x,t,K,K1,K2,K3,K4;

const double sixth = 1.0 / 6.0;

static double dx_dt(double t, double x){
return tan(x) + 1;
}

int main(int argc, const char * argv[]) {

/*======================================================================*/
/*===================== Runge-Kutta Method for ODE =====================*/
/*======================================================================*/

double t_initial = 1.0;// initial time
double x_initial = 1.0;// initial x position

double t_final = 1.1;// value of t wish to know x
double dt = 0.025;// time interval for updates
double halfdt = 0.5*dt;

/*======================================================================*/

while(t_initial < t_final){

/*============================ Runge-Kutta increments =================================*/

double K1 = dt*dx_dt( t_initial, x_initial );
double K2 = dt*dx_dt( t_initial + halfdt, x_initial + halfdt*K1 );
double K3 = dt*dx_dt( t_initial + halfdt, x_initial + halfdt*K2 );
double K4 = dt*dx_dt( t_initial + dt, x_initial + dt*K3 );

x_initial += sixth*(K1 + 2*(K2 + K3) + K4);

/*============================ prints =================================*/

std::cout << t_initial << std::setw(16) << x_initial << "\n";

/*============================ re-setting update conditions =================================*/

t_initial += dt;

/*======================================================================*/
}

std::cout<<"----------------------------------------------\n";
std::cout << "t =  "<< t_initial << ",  x =  "<< x_initial << std::endl;}/* main */

4

c++math numerical-methods runge-kutta

Решение

Другие решения

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

Должно быть либо

k2 = dt*f(t+0.5*dt, x+0.5*k1)

или же

k2 = f(t+0.5*dt, x+0.5*dt*k1)

другой kс соответственно.

Обратите внимание, что в обоих случаях наклон fтолько умножается на dt один раз.

4

Я думаю, что вы включили слишком много приращений и создали проблемы, переставив математику. Попробуй это:

#include <math.h>
#include <iostream>
#include <iomanip>

static double dx_dt(double t, double x)
{
return tan(x) + 1;
}

int main(int argc, const char * argv[])
{
double t = 1.0;
double t_end = 1.1;

double y = 1.0;
double h = 0.025;

std::cout << std::setprecision(16);

int n = static_cast<int>((t_end - t) / h);

for (int i = 0; i < n; i++)
{
double k1 = dx_dt(t, y);
double k2 = dx_dt(t + h / 2.0, y + h*k1 / 2.0);
double k3 = dx_dt(t + h / 2.0, y + h*k2 / 2.0);
double k4 = dx_dt(t + h, y + h*k3);

y += (k1 + 2 * k2 + 2 * k3 + k4) * h / 6.0;

std::cout << t << ": " << y << std::endl;

t += h;
}

std::cout << "----------------------------------------------\n";
std::cout << "t =  " << t << ",  x =  " << y << std::endl;

std::getchar();
}

Я предварительно вычисляю, сколько раз выполнить итерацию, это позволяет избежать нескольких разных проблем. Также, как уже упоминали другие, работающий пример в Википедии для двухэтапного варианта алгоритма.

Я взял на себя смелость изменить имена переменных в соответствии с Википедией. Хороший совет всегда соответствует названию вашего ссылочного текста, пока вещь не сработает.

3

Источник

Accepted Answer

Проблема в том, что таблица, используемая для вашего кода, отличается от таблицы, которую вы разместили в википедии. Тот, который вы используете, это:

0
1/2     1/2
1/2     0       1/2
1       0       0       1
1/6     1/3     1/3     1/6

И тот, который используется в Википедии

0
2/3     2/3
1/4     3/4

Различные таблицы будут давать разные результаты в зависимости от размера шага, который используется для убедитесь, что размер шага достаточно хорош для определенной точности. Однако когда dt -> 0тогда все столы одинаковы.

Помимо всего этого, ваш код в любом случае неверен даже для RK4. Вторая часть функции должна иметь половинки, а не 0.5*dt:

double K1 = dt*dx_dt( t_initial, x_initial );
double K2 = dt*dx_dt( t_initial + halfdt, x_initial + 0.5*K1 );
double K3 = dt*dx_dt( t_initial + halfdt, x_initial + 0.5*K2 );
double K4 = dt*dx_dt( t_initial + dt, x_initial + K3 );

4