Рунге-Кутт 4 порядка, диффур 2-го порядка

[sko1111]

\[y'' = -2\left(\frac{x-a}{\sigma}\right)y' - \frac 2\sigma y\]

\[\begin{split} y(0) = 0,1\\ y'(0) = -2\\ a = 2\\ \sigma = 0,2 \end{split}\]

преобразовал в систему:

\[\left\{\begin{split} U' = f(x, U, y) = -2\left(\frac{x-a}{\sigma}\right)U - \frac 2\sigma y\\ y' = g(x, U, y) = U \end{split}\right.\]

реализовал на C:

#include <stdio.h>

double a = 2.0;
double sigma = 0.2;

double f(double x, double u, double y) {
   return -2*u*(x-a)/sigma-2*y/sigma;
}

double g(double x, double u, double y) {
   return u;
}

int main() {
   double yn = 0.1;
   double ysn = -2.0;
   double x = 0.0f;
   double h = 1e-3;
   double b = 0.5;

   double k1,k2,k3,k4,q1,q2,q3,q4;

   while (x <= b+h) {
      k1 = f(x, ysn, yn)*h;
      q1 = g(x, ysn, yn)*h;
      k2 = f(x+h/2, ysn+k1/2, yn+q1/2)*h;
      q2 = g(x+h/2, ysn+k1/2, yn+q1/2)*h;
      k3 = f(x+h/2, ysn+k2/2, yn+q2/2)*h;
      q3 = g(x+h/2, ysn+k2/2, yn+q2/2)*h;
      k4 = f(x+h, ysn+k3, yn+q3)*h;
      q4 = g(x+h, ysn+k3, yn+q3)*h;
      ysn += (k1+2*k2+2*k3+k4)/6;
      yn += (q1+2*q2+2*q3+q4)/6;
      printf("%e\n",(k1+2*k2+2*k3+k4)/6-f(x, ysn, yn));
      x += h;
   }
   return 0;
}

printf по идее должен выводить невязку:

Цитата:

4.176625e+01
4.256773e+01
4.338410e+01
4.421562e+01
4.506257e+01
4.592522e+01
4.680385e+01
4.769874e+01
4.861020e+01
4.953852e+01
5.048399e+01
5.144692e+01
5.242764e+01
5.342644e+01

........................

1.317720e+05
1.335851e+05
1.354217e+05
1.372820e+05
1.391664e+05
1.410750e+05
1.430083e+05
1.449664e+05

подскажите, где неверно? я уже не знаю где искать

[Lotos]

А что Вы выводите на экран? Я не специалист по численным методам, но метод Р.-К. использовала не один раз?
Может в выводе при вызове функции f надо аргумент x+h?

[sko1111]

ну уравнение имеет вид

\[y'' = f(x,y,y')\]

. я вычисляю невязку подстановкой численного решения в уравнение.
x+h --- те же яйца, вид сбоку тоже e+5 выходит.

как это так у всех получается, а у меня нет

[V.V.]

У Вас все $k_i$ умножены на $h$. Поэтому фактически Вы выводите $-f$.

[sko1111]

вообще там схема такая (по моим лекциям):

\[\begin{split} k_1 = f(x, y, z)\\ q_1 = g(x, y, z)\\ k_2 = f(x+\frac h2, y+\frac{hk_1}2, z+\frac{hq_1}2)\\ q_2 = g(x+\frac h2, y+\frac{hk_1}2, z+\frac{hq_1}2)\\ k_3 = f(x+\frac h2, y+\frac{hk_2}2, z+\frac{hq_2}2)\\ q_3 = g(x+\frac h2, y+\frac{hk_2}2, z+\frac{hq_2}2)\\ k_4 = f(x+h, y+hk_3, z+hq_3)\\ q_4 = g(x+h, y+hk_3, z+hq_3)\\ \end{split}\]

просто я внес $h$ в коэффициенты. в гугле смотрел, там так же.
а не выходит.
можно поподробнее, почему я вывожу $-f$? ведь при выводе я на $h$ уже не домножаю. функция слишком быстро растет, не похоже на верный ответ.

[Lotos]

А вы пробовали выводить при f(x+h,....)? Уравнение зависит от х напрямую?
Получается е+01
А еще надо просто тестировать программу на уравнении, решение которого вы можете аналитически определить. Решаете уравнение аналитически и находите значение y(x) с помощью программы, это же уравнение решаете численно и сравниваете результаты.
Правда, есть еще нюансы с устойчивостью метода для вашего уравнения, но это уже тонкости...
Откуда такая уверенность, что программа написана неправильно?