У меня есть следующий алгоритм 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);
}
почему-то... -- Что такое отладчик?. Вам нужно выполнить расчет вручную, а затем сравнить результаты при отладке вашей программы, чтобы увидеть, где (и когда) ваша реализация отличается от того, что вы рассчитали вручную.
В вашей реализации Рунге-Кутты нет ничего плохого — просто ваше ожидание решения.
Немного переформулируйте уравнение и используйте интегрирующий коэффициент. Вы найдете точное решение
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
Вы можете проверить простым дифференцированием, что это решает dy/dx=x^2+y и проходит через точку (x0,y0). Ваша константа c1 равна y0+x0^2+2x0+2. Это явно зависит от x0. В противном случае вы бы не достигли граничного условия.
Если вы установите x0=-5, то вы применяете начальное граничное условие y(-5)=-2, а НЕ y(0)=-2. Так что ваш сценарий WolframAlpha не имеет значения.
Оооооооооооооооооооооооооооооо ты так прав
какой у Вас вопрос?