Численное решение волнового уравнения

[lv00]

Здравствуйте. Я пытаюсь сделать простую задачу - решить однородное линейное волновое уравнение

$$u_{tt} = v^2 u_{xx}$$

,

где $v$ - скорость распространяющейся волны.
В качестве начальных условий я использую

$$u(x, 0) = \sin(x)$$

$$u_t(x, 0) = -v \sin(x) $$

что соответствует волне с профилем $\sin(x-vt)$.

В течение первых двух-трех итераций полученное решение с помощью метода Эйлера (Рунге-Кутты 2го порядка точности) совпадает с аналитическим ($\sin(x-vt)$), но далее происходит резкий "выстрел" - значение решения уходит на отметку около 128. В чем причина - я не могу понять.

Привожу код, написанный на C++.

//2nd derivative function

    float z2der(int n, int i, int j, float* x, float** zx) { 
   float h = x[2] - x[1];
   float zr, zl, zm, zrr, zll;

   zm = -30.0*zx[i][j];
   if (i<n - 1)
      zr = 16.0*zx[i + 1][j];
   else
      zr = 16.0;
   if (i>0)
      zl = 16.0*zx[i - 1][j];
   if (i<n - 2)
      zrr = -zx[i + 2][j];
   if (i>1)
      zll = -zx[i - 2][j];
   
   if (i>1 && i<n-2)
      return (zll + zrr + zl + zr + zm) / (12.0 * h * h);
   if (i == 1)
      return (15 / 4 * zx[i][j] - 77 / 6 * zx[i + 1][j] - 107 / 6 * zx[i + 2][j] - 13 * zx[i + 3][j] - 61 / 12 * zx[i + 4][j]) / (h*h);
   if (i==n-2)
      return (15 / 4 * zx[i][j] - 77 / 6 * zx[i - 1][j] - 107 / 6 * zx[i - 2][j] - 13 * zx[i - 3][j] - 61 / 12 * zx[i - 4][j]) / (h*h);
   if (i == 0 || i == n - 1)
      return 0;
}


//PARAMETERS

    float tmin = 0.0;
   float tmax = 12.0;
   int nt = 37000;
   float dt = (tmax - tmin) / float(nt);
   float *time = new float[nt];
   for (int j = 0; j<nt; j++)
      time[j] = tmin + dt*float(j);

   float xmin = -30.0;
   float xmax = 30.0;
   int nx = 901;
   float *x = new float[nx];
   for (int i = 0; i<nx; i++)
      x[i] = xmin + (xmax - xmin)*float(i) / float(nx);

   float v = -0.1; // velocity

    float** wave = new float*[nx];
   for (int i = 0; i < nx; i++)
      wave[i] = new float[nt];

   float** waved = new float*[nx];
   for (int i = 0; i < nx; i++)
      waved[i] = new float[nt];

   for (int i = 0; i<nx; i++) {
      xx = x[i];
      wave[i][0] = sin(xx);
      waved[i][0] = -v*v*sin(xx);
   }

   float* wavedd = new float[nx];
   float* waveddot = new float[nx];

   for (int i = 0; i<nx; i++) { //j=0;
      wavedd[i] = z2der(nx, i, 0, x, wave);
      waveddot[i] = v*v*wavedd[i];
   }

   for (int j = 0; j<nt - 1; j++) {
      for (int i = 0; i<nx; i++) {
         waved[i][j + 1] = waved[i][j] + dt * waveddot[i];
         wave[i][j + 1] = wave[i][j] + dt * waved[i][j];

         wavedd[i] = z2der(nx, i, j + 1, x, wave);
         waveddot[i] = v*v*wavedd[i];
      }
   }

[Pphantom]

Ну, помимо того, что код в целом ужасен... Как Вы думаете, чему равно, например, 15 / 4 * zx[i][j]?

(Ответ)

А вот и не угадали, 3 * zx[i][j].

[Red_Herring]

Кстати, явная схема для волнового уравнения условно устойчива

[Pavia]

Red_Herring

Red_Herring в сообщении #1044247 писал(а):

Кстати, явная схема для волнового уравнения условно устойчива

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


PS. На самом деле у меня тоже затык с численным решением. Но пока ещё сам по разбираюсь.

[Red_Herring]

Pavia в сообщении #1044344 писал(а):

Что вы имеете в виду?

Ну то что если шаги по времени и пространству $\tau,h$ то при $v\tau >h$ устойчивость в явной (естественной) схеме теряется.

[Pavia]

Red_Herring
Обычно даётся условие устойчивости в общем виде. И оно работает для любой разностной схемы (явной, не явной).
$\frac{v^2\tau^2}{h^2} <1$
Если у нас явная схема, то $h>0$, $\tau>0$, но упростить всё равно не получится из-за скорости. Хотя если её взять по модулю, то можно.

Red_Herring в сообщении #1044247 писал(а):

Кстати, явная схема для волнового уравнения условно устойчива

У автора в коде есть следующие параметры.

   float tmin = 0.0;
   float tmax = 12.0;
   int nt = 37000;
   float dt = (tmax - tmin) / float(nt);

   float xmin = -30.0;
   float xmax = 30.0;
   int nx = 901;

 float v = -0.1; // velocity

Откуда имеем, данные:
$\tau=3,24 \cdot 10^{-4}$
$h=0,066$
$v=-0.1$

$\frac{0.01 \cdot 1,049 \cdot 10^{-7}}{0.0044}=2.385 \cdot 10^{-8}<1$
Условие верно.

[lv00]

Pphantom в сообщении #1044137 писал(а):

Ну, помимо того, что код в целом ужасен... Как Вы думаете, чему равно, например, 15 / 4 * zx[i][j]?

(Ответ)

А вот и не угадали, 3 * zx[i][j].

Данная ошибка у меня сделана только в определении производной на границе. Сейчас я ее исправил, но все равно результат аналогичный - на третьей итерации решение "взлетает".

P.S. В чем ужасность кода?

[Pavia]

lv00 в сообщении #1044362 писал(а):

P.S. В чем ужасность кода?

Ваш код понять стороннему человеку просто невозможно.
Разве из вашего кода можно понять вы используете метод Эйлера или метод Рунге-Кутты?
Много к чему можно придраться
1) Функций мало.
2) Комментариев мало.
3) Имена переменных непонятны.
4) Большинство констант не понятно откуда взялись.

zm = -30.0*zx[i][j];

Откуда тут -30?

  zm = -30.0*zx[i][j];
   if (i<n - 1)
      zr = 16.0*zx[i + 1][j]; // ';' - лишняя
   else
      zr = 16.0;

Много if в первой функции. За ними теряется суть. Их надо вынести из цикла.
Что-бы не писать много if для обработки границ,
надо просто увеличив массив с каждой стороны. И заполнить приближенными значениями, к примеру просто продублировать крайние значения.

[Pphantom]

Red_Herring в сообщении #1044346 писал(а):

Ну то что если шаги по времени и пространству $\tau,h$ то при $v\tau >h$ устойчивость в явной (естественной) схеме теряется.

Тут есть конкретные параметры, для них это условие выполняется. Для распада на 2-3 шаге это не причина.

Больше похоже на проявление немонотонности, но ввиду жуткого и при этом неполного кода это и глазами проверить трудно, и запустить и посмотреть нельзя.

lv00 в сообщении #1044362 писал(а):

P.S. В чем ужасность кода?

Из него совершенно непонятно, что и как Вы делаете. Комментариев нет, расшифровки названий переменных нет, везде напиханы "магические константы"...

[arseniiv]

Pavia в сообщении #1044397 писал(а):

// ';' - лишняя

Ничего не лишняя. Она входит в оператор выражение ; и не входит в оператор if ( выражение ) оператор else оператор.

[NT2000]

lv00 в сообщении #1044362 писал(а):

на третьей итерации решение "взлетает".

Я тут быстро поправил огрехи и скомпилировал ваш код, получил данные
для первых 5 итераций:

i=0; sinus(-30)=0.988032; -v*v*sin(xx)=-0.00988032
i=1; sinus(-29.9334)=0.996106; -v*v*sin(xx)=-0.00996106
i=2; sinus(-29.8668)=0.999765; -v*v*sin(xx)=-0.00999765
i=3; sinus(-29.8002)=0.998992; -v*v*sin(xx)=-0.00998992
i=4; sinus(-29.7336)=0.99379; -v*v*sin(xx)=-0.0099379
i=5; sinus(-29.667)=0.984183; -v*v*sin(xx)=-0.00984183

i=0; wavedd[0] = 0
i=1; wavedd[1] = -10110.2
i=2; wavedd[2] = -0.999792
i=3; wavedd[3] = -0.999039
i=4; wavedd[4] = -0.993842
i=5; wavedd[5] = -0.984128

i=0, j=0; waveddot[0] = 0
i=1, j=0; waveddot[1] = 8.42366
i=2, j=0; waveddot[2] = -2.827
i=3, j=0; waveddot[3] = -2.81316
i=4, j=0; waveddot[4] = -2.78684
i=5, j=0; waveddot[5] = -2.74817
i=0, j=1; waveddot[0] = 0
i=1, j=1; waveddot[1] = 8.42355
i=2, j=1; waveddot[2] = -2.82702
i=3, j=1; waveddot[3] = -2.81314
i=4, j=1; waveddot[4] = -2.78683
i=5, j=1; waveddot[5] = -2.74817
i=0, j=2; waveddot[0] = 0
i=1, j=2; waveddot[1] = 8.42344
i=2, j=2; waveddot[2] = -2.82704
i=3, j=2; waveddot[3] = -2.81313
i=4, j=2; waveddot[4] = -2.78682
i=5, j=2; waveddot[5] = -2.74816
i=0, j=3; waveddot[0] = 0
i=1, j=3; waveddot[1] = 8.42333
i=2, j=3; waveddot[2] = -2.82705
i=3, j=3; waveddot[3] = -2.81312
i=4, j=3; waveddot[4] = -2.78681
i=5, j=3; waveddot[5] = -2.74815

Скажите о каких взлетах вы говорите?

[lv00]

Привожу исправленный код. If-ы не изменял, добавил комментарии под переменные. Константы - такие как 30.0, 16.0 и пр. - это коэффициенты для аппроксимации второй производной конечной разностью 4го порядка точности.

    float 2der(int n, int i, int j, float* x, float** func) { //Second Der. of func
   float h = x[2] - x[1];
   float zr, zl, zm, zrr, zll;

   zm = -30.0*func[i][j];
   if (i<n - 1)
      zr = 16.0*func[i + 1][j];
   if (i>0)
      zl = 16.0*func[i - 1][j];
   if (i<n - 2)
      zrr = -func[i + 2][j];
   if (i>1)
      zll = -func[i - 2][j];
   
   if (i>1 && i<n-2)
      return (zll + zrr + zl + zr + zm) / (12.0 * h * h); 
   if (i == 1)
      return (15.0 / 4.0 * func[i][j] - 77.0 / 6.0 * func[i + 1][j] - 107.0 / 6.0 * func[i + 2][j] - 13.0 * func[i + 3][j] - 61.0 / 12.0 * func[i + 4][j]) / (h*h);
   if (i==n-2)
      return (15.0 / 4.0 * func[i][j] - 77.0 / 6.0 * func[i - 1][j] - 107.0 / 6.0 * func[i - 2][j] - 13.0 * func[i - 3][j] - 61.0 / 12.0 * func[i - 4][j]) / (h*h);
   if (i == 0 || i == n - 1)
      return 0;
}


//PARAMETERS

    float tmin = 0.0;
   float tmax = 12.0;
   int nt = 37000; //number of points for time
   float dt = (tmax - tmin) / float(nt);
   float *time = new float[nt];
   for (int j = 0; j<nt; j++)
      time[j] = tmin + dt*float(j);

   float xmin = -30.0;
   float xmax = 30.0;
   int nx = 901; //number of points for x-coordinate
   float *x = new float[nx];
   for (int i = 0; i<nx; i++)
      x[i] = xmin + (xmax - xmin)*float(i) / float(nx);

   float v = -0.1; // velocity of wave

//SOLUTION

    float** wave = new float*[nx]; //array of solution
   for (int i = 0; i < nx; i++)
      wave[i] = new float[nt];

   float** waved = new float*[nx]; //array of 1st derivative of solution
   for (int i = 0; i < nx; i++)
      waved[i] = new float[nt];

   for (int i = 0; i<nx; i++) { //initial conditions
      xx = x[i];
      wave[i][0] = sin(xx);
      waved[i][0] = -v*v*sin(xx);
   }

   float* wavedd = new float[nx]; //auxiliary variables
   float* waveddot = new float[nx];

   for (int i = 0; i<nx; i++) { //j=0;
      wavedd[i] = 2der(nx, i, 0, x, wave);
      waveddot[i] = v*v*wavedd[i];
   }

   for (int j = 0; j<nt - 1; j++) { // Euler method solution
      for (int i = 0; i<nx; i++) {
         waved[i][j + 1] = waved[i][j] + dt * waveddot[i];
         wave[i][j + 1] = wave[i][j] + dt * waved[i][j];

         wavedd[i] = 2der(nx, i, j + 1, x, wave);
         waveddot[i] = v*v*wavedd[i];
      }
   }

Pavia в сообщении #1044397 писал(а):

lv00 в сообщении #1044362 писал(а):

P.S. В чем ужасность кода?

Откуда тут -30?

-30 тут, т.к. это коэффициент для разностной схемы, аппроксимирующий вторую производную. Там должно быть -$5/2$, как раз $-30/12 = -5/2$.

NT2000 в сообщении #1044455 писал(а):

lv00 в сообщении #1044362 писал(а):

на третьей итерации решение "взлетает".

Скажите о каких взлетах вы говорите?

Вы написали итерации по координате (соответствует i). Я же говорю про итерации по времени (как раз работа алгоритма численного решения Рунге-Кутты), они соответствуют j.

[Pavia]

lv00 в сообщении #1044677 писал(а):

if (i == 1)
return (15.0 / 4.0 * func[i][j] - 77.0 / 6.0 * func[i + 1][j] - 107.0 / 6.0 * func[i + 2][j] - 13.0 * func[i + 3][j] - 61.0 / 12.0 * func[i + 4][j]) / (h*h);
if (i==n-2)
return (15.0 / 4.0 * func[i][j] - 77.0 / 6.0 * func[i - 1][j] - 107.0 / 6.0 * func[i - 2][j] - 13.0 * func[i - 3][j] - 61.0 / 12.0 * func[i - 4][j]) / (h*h);

Вот тут ошибка. Если проверить по норме.
$5.0 / 4.0 * 1 - 77.0 / 6.0 * 1 - 107.0 / 6.0 * 1 - 13.0 * 1 - 61.0 / 12.0 *1=-47,5$
А должен быть 0.

Цитата:

if (i == 0 || i == n - 1)
return 0;

Не лучшее решение.

[lv00]

Pavia в сообщении #1044750 писал(а):

lv00 в сообщении #1044677 писал(а):

if (i == 1)
return (15.0 / 4.0 * func[i][j] - 77.0 / 6.0 * func[i + 1][j] - 107.0 / 6.0 * func[i + 2][j] - 13.0 * func[i + 3][j] - 61.0 / 12.0 * func[i + 4][j]) / (h*h);
if (i==n-2)
return (15.0 / 4.0 * func[i][j] - 77.0 / 6.0 * func[i - 1][j] - 107.0 / 6.0 * func[i - 2][j] - 13.0 * func[i - 3][j] - 61.0 / 12.0 * func[i - 4][j]) / (h*h);

Вот тут ошибка. Если проверить по норме.
$5.0 / 4.0 * 1 - 77.0 / 6.0 * 1 - 107.0 / 6.0 * 1 - 13.0 * 1 - 61.0 / 12.0 *1=-47,5$
А должен быть 0.

Я поправил. Правильные коэффициенты: $15.0/4.0, -77.0/6.0, 107.0/6.0, -13.0, 61.0/12.0, -5.0/6.0$

Цитата:

Цитата:

if (i == 0 || i == n - 1)
return 0;

Не лучшее решение.

Почему?

После исправления все равно наблюдаются "взлеты". Привожу данные.

i=0: wave[0][0] = 0.988032, wave[0][1] = 0.988028, wave[0][2] = 0.988025, wave[0][3] = 0.988022
i=1: wave[1][0] = 0.996106, wave[1][1] = 0.996103, wave[1][2] = 0.9961, wave[1][3] = 384.915
i=2: wave[2][0] = 0.999765, wave[2][1] = 0.999762, wave[2][2] = 0.999758, wave[2][3] = -126.973
i=3: wave[3][0] = 0.998992, wave[3][1] = 0.998989, wave[3][2] = 0.998985, wave[3][3] = -126.974

и так далее.

[lv00]

UPD: заметил еще опечатку в начальном условии для производной по времени - должно быть $v \cos(xx)$, а не $-v^2 \sin(xx)$, но это не влияет на ситуацию - ровно такой же рост в том же месте.