2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




На страницу Пред.  1, 2
 
 Re: Численное решение волнового уравнения
Сообщение12.08.2015, 15:59 
Аватара пользователя
Внес ваши изменения.
Ну никак не пойму, где фиксировать указанный взлет?
Ткните пальцем.
Код:
i=0; sinus(-30)=0.988032; v*cos(xx)=-0.0154251
i=1; sinus(-29.9334)=0.996106; v*cos(xx)=-0.0088162
i=2; sinus(-29.8668)=0.999765; v*cos(xx)=-0.00216837
i=3; sinus(-29.8002)=0.998992; v*cos(xx)=0.00448927
i=4; sinus(-29.7336)=0.99379; v*cos(xx)=0.011127
i=5; sinus(-29.667)=0.984183; v*cos(xx)=0.0177154

i=0; wavedd[0] = 0
i=1; wavedd[1] = -10292.5
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; wave[0,0] = 0.988032
i=1, j=0; wave[1,0] = 0.996106
i=2, j=0; wave[2,0] = 0.999765
i=3, j=0; wave[3,0] = 0.998992
i=4, j=0; wave[4,0] = 0.99379
i=5, j=0; wave[5,0] = 0.984183
i=0, j=1; wave[0,1] = 0.988027
i=1, j=1; wave[1,1] = 0.996103
i=2, j=1; wave[2,1] = 0.999764
i=3, j=1; wave[3,1] = 0.998993
i=4, j=1; wave[4,1] = 0.993794
i=5, j=1; wave[5,1] = 0.984189
i=0, j=2; wave[0,2] = 0.988022
i=1, j=2; wave[1,2] = 0.996209
i=2, j=2; wave[2,2] = 0.999763
i=3, j=2; wave[3,2] = 0.998995
i=4, j=2; wave[4,2] = 0.993797
i=5, j=2; wave[5,2] = 0.984195
i=0, j=3; wave[0,3] = 0.988017
i=1, j=3; wave[1,3] = 0.996305
i=2, j=3; wave[2,3] = 0.999766
i=3, j=3; wave[3,3] = 0.998999
i=4, j=3; wave[4,3] = 0.993804
i=5, j=3; wave[5,3] = 0.984203
i=0, j=4; wave[0,4] = 0.988012
i=1, j=4; wave[1,4] = 0.996393
i=2, j=4; wave[2,4] = 0.999771
i=3, j=4; wave[3,4] = 0.999007
i=4, j=4; wave[4,4] = 0.993814
i=5, j=4; wave[5,4] = 0.984215
i=0, j=5; wave[0,5] = 0.988007
i=1, j=5; wave[1,5] = 0.996472
i=2, j=5; wave[2,5] = 0.999779
i=3, j=5; wave[3,5] = 0.999017
i=4, j=5; wave[4,5] = 0.993826
i=5, j=5; wave[5,5] = 0.984229


PS. Вот еще замечание.
Код:
   for (int i = 0; i<nx; i++) { //initial conditions
      xx = x[i];              //Вот первое использование переменной ХХ,
                              //а где она определена, какой имеет тип?

Пишу об этом, потому как 2-ой раз приводите код, а определения переменной как не было, так и нет.

 
 
 
 Re: Численное решение волнового уравнения
Сообщение12.08.2015, 17:15 
Да, у вас все нормально. Но у меня почему-то, например, на
Код:
i=1, j=3: wave[1][3] = 384.915

что меня и беспокоит (все результаты до шага j=3 совпадают с вашими).
Переменная xx определена следующим образом:
Код:
float xx;

 
 
 
 Re: Численное решение волнового уравнения
Сообщение12.08.2015, 17:46 
Аватара пользователя
lv00 в сообщении #1044817 писал(а):
Да, у вас все нормально. Но у меня почему-то, например, на
Код:
i=1, j=3: wave[1][3] = 384.915

что меня и беспокоит (все результаты до шага j=3 совпадают с вашими).

Мне не остается как опубликовать код который использую.

(Оффтоп)

Код:
#include <cstdlib>
#include <iostream>
#include <math.h>

using namespace std;

//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.0/4.0 * zx[i][j] - 77.0/6.0 * zx[i + 1][j] - 107.0/6.0 * zx[i + 2][j] - 13.0 * zx[i + 3][j] - 61.0/12.0 * zx[i + 4][j] - 5.0/6.0 * zx[i + 5][j]) / (h*h);
   if (i==n-2)
      return (15.0/4.0 * zx[i][j] - 77.0/6.0 * zx[i - 1][j] - 107.0/6.0 * zx[i - 2][j] - 13.0 * zx[i - 3][j] - 61.0/12.0 * zx[i - 4][j]) - 5.0/6.0 * zx[i - 5][j]/ (h*h);
   if (i == 0 || i == n - 1)
      return 0;
}

int main(){
//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
      float xx = x[i];
      wave[i][0] = sin(xx);
      waved[i][0] = v*cos(xx);
      if (i<6){
         cout<< "i="<< i << "; sinus(" << xx << ")="<< sin(xx) << "; ";
         cout<< "v*cos(xx)=" <<  v*cos(xx) << endl;
      }
   }

   float* wavedd = new float[nx];   //auxiliary variables
   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*wavedd[i];
    if (i<6){
         cout<< "i="<< i << "; wavedd[" << i <<"] = "<< wavedd[i]  << endl;
      }
   }

   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];
         
      if ((i<6)&&(j<6)){
            cout<< "i="<< i <<", j=" << j << "; wave[" << i <<"," << j <<"] = "<<wave[i][j]<< endl;
         }

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

 
 
 
 Re: Численное решение волнового уравнения
Сообщение12.08.2015, 20:39 

(Оффтоп)

Друзья, давайте жить дружно! Давайте использовать тег [​syntax]! Для C++ подсветка синтаксиса есть. (Если нужны подробности, см. здесь.)

P. S. Заодно там и код сворачивается автоматически.

 
 
 
 Re: Численное решение волнового уравнения
Сообщение12.08.2015, 21:08 
Аватара пользователя
( symntax :D )
Виноват. Я такое искал в теге [code], теперь буду знать.
Повторим код, для закрепления.
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
#include <cstdlib>
#include <iostream>
#include <math.h>

using namespace std;

//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.0/4.0 * zx[i][j] - 77.0/6.0 * zx[i + 1][j] - 107.0/6.0 * zx[i + 2][j] - 13.0 * zx[i + 3][j] - 61.0/12.0 * zx[i + 4][j] - 5.0/6.0 * zx[i + 5][j]) / (h*h);
   if (i==n-2)
      return (15.0/4.0 * zx[i][j] - 77.0/6.0 * zx[i - 1][j] - 107.0/6.0 * zx[i - 2][j] - 13.0 * zx[i - 3][j] - 61.0/12.0 * zx[i - 4][j]) - 5.0/6.0 * zx[i - 5][j]/ (h*h);
   if (i == 0 || i == n - 1)
      return 0;
}

int main(){
//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
      float xx = x[i];
      wave[i][0] = sin(xx);
      waved[i][0] = v*cos(xx);
      if (i<6){
         cout<< "i="<< i << "; sin(" << xx << ")="<< sin(xx) << "; ";
         cout<< "v*cos(xx)=" <<  v*cos(xx) << endl;
      }
   }

   float* wavedd = new float[nx];   //auxiliary variables
   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*wavedd[i];
    if (i<6){
         cout<< "i="<< i << "; wavedd[" << i <<"] = "<< wavedd[i]  << endl;
      }
   }

   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];
         
         if ((i<6)&&(j<6)){
            cout<< "i="<< i <<", j=" << j << "; wave[" << i <<"," << j <<"] = "<<wave[i][j]<< endl;
         }

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

 
 
 
 Re: Численное решение волнового уравнения
Сообщение12.08.2015, 21:54 

(Ой.)

NT2000 в сообщении #1044850 писал(а):
( symntax :D )
Поправил, спасибо!

 
 
 
 Re: Численное решение волнового уравнения
Сообщение13.08.2015, 21:25 
NT2000 в сообщении #1044850 писал(а):
( symntax :D )
Виноват. Я такое искал в теге code, теперь буду знать.
Повторим код, для закрепления.
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
#include <cstdlib>
#include <iostream>
#include <math.h>

using namespace std;

//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.0/4.0 * zx[i][j] - 77.0/6.0 * zx[i + 1][j] - 107.0/6.0 * zx[i + 2][j] - 13.0 * zx[i + 3][j] - 61.0/12.0 * zx[i + 4][j] - 5.0/6.0 * zx[i + 5][j]) / (h*h);
   if (i==n-2)
      return (15.0/4.0 * zx[i][j] - 77.0/6.0 * zx[i - 1][j] - 107.0/6.0 * zx[i - 2][j] - 13.0 * zx[i - 3][j] - 61.0/12.0 * zx[i - 4][j]) - 5.0/6.0 * zx[i - 5][j]/ (h*h);
   if (i == 0 || i == n - 1)
      return 0;
}

int main(){
//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
      float xx = x[i];
      wave[i][0] = sin(xx);
      waved[i][0] = v*cos(xx);
      if (i<6){
         cout<< "i="<< i << "; sin(" << xx << ")="<< sin(xx) << "; ";
         cout<< "v*cos(xx)=" <<  v*cos(xx) << endl;
      }
   }

   float* wavedd = new float[nx];   //auxiliary variables
   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*wavedd[i];
    if (i<6){
         cout<< "i="<< i << "; wavedd[" << i <<"] = "<< wavedd[i]  << endl;
      }
   }

   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];
         
         if ((i<6)&&(j<6)){
            cout<< "i="<< i <<", j=" << j << "; wave[" << i <<"," << j <<"] = "<<wave[i][j]<< endl;
         }

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


Я скопировал ваш код, и вот что получил:
например
Код:
wave[1,1] = 0.996103
wave[1,2] = 0.996209
wave[1,3] = -50761.6

и дальше все больше и больше с каждой итерацией по времени.

Компилирую в Visual Studio 2013.

 
 
 
 Re: Численное решение волнового уравнения
Сообщение13.08.2015, 22:57 
Аватара пользователя
lv00
Есть такая вещь как отладчик. Возьмите код и пройдитесь по шагам.
Всё же легко вычисляется.

Код:
        wavedd[i] = z2der(nx, i, j + 1, x, wave);

Откуда тут +1 взялось? Остальные ошибки ищите сами.

 
 
 
 Re: Численное решение волнового уравнения
Сообщение14.08.2015, 00:46 
Аватара пользователя
Pavia в сообщении #1045148 писал(а):
Откуда тут +1 взялось? Остальные ошибки ищите сами.

Да - всё верно :appl:
Повторил компиляцию под VS2012 с вызовом функции :
Код:
        wavedd[i] = z2der(nx, i, j + 1, x, wave);  //фиксируется взлет

i=1, j=3; wave[1,3] = -50761.6
i=2, j=3; wave[2,3] = 1280.73
i=3, j=3; wave[3,3] = 1280.73
i=4, j=3; wave[4,3] = 1280.72
i=5, j=3; wave[5,3] = 1280.71
i=0, j=4; wave[0,4] = 0.988012
i=1, j=4; wave[1,4] = -152287
i=2, j=4; wave[2,4] = 3840.19
i=3, j=4; wave[3,4] = 3840.19
i=4, j=4; wave[4,4] = 3840.18

Если поправить - то ок.

PS.Что интересно, для компилятора Dev-C++ - в обоих случаях есть ок.

 
 
 
 Re: Численное решение волнового уравнения
Сообщение14.08.2015, 15:29 
Аватара пользователя
Мдаа...
Под VS2012 - проверил и иногда программа сбоила даже с правильным вызовом функции z2der.
Мне кажется всё дело в приделении памяти через операторы new float[].
Динамическая аллокация памяти - только выделяет непрерывную область
для переменной указанного типа, но не инициилирует её нулями.
А ТС пробует эту область с мусором использовать в вычислениях imho.

 
 
 
 Re: Численное решение волнового уравнения
Сообщение14.08.2015, 16:13 

(Ещё немного о тегах.)

А ещё есть [​tt], оформляющий код и цветом (зелёный, как и у [​code]), и моноширинностью: ABC123. Есть кнопка над полем набора, как и для остальных. :-)

 
 
 
 Re: Численное решение волнового уравнения
Сообщение15.08.2015, 20:24 
Аватара пользователя
Думаю ТС уже всё поправил, коль скоро ничего не пишет,
но всё равно опубликую свою версию - одинаково работающую независимо от компилятора.
Добавлено инициация выделенной области, плюс поправки в знаках в граничных точках.
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
  1. // ConsoleApplication1.cpp
  2. //
  3. #include <iostream>
  4. #include <math.h>
  5. #include "stdafx.h"
  6.  
  7. using namespace std;
  8.  
  9. //2nd derivative function
  10. float z2der(int n, int i, int j, float* x, float** zx) {
  11.    float h = x[2] - x[1];
  12.    float zr, zl, zm, zrr, zll;
  13.  
  14.    zm = -30.0f*zx[i][j];
  15.    if (i<n - 1)
  16.       zr = 16.0f*zx[i + 1][j];
  17.    else
  18.       zr = 16.0f;
  19.    if (i>0)
  20.       zl = 16.0f*zx[i - 1][j];
  21.    if (i<n - 2)
  22.       zrr = -zx[i + 2][j];
  23.    if (i>1)
  24.       zll = -zx[i - 2][j];
  25.    
  26.    if (i>1 && i<n-2)
  27.       return (zll + zrr + zl + zr + zm) / (12.0f * h * h);
  28.    if (i == 1)
  29.       return (15.0f/4.0f * zx[i][j] - 77.0f/6.0f * zx[i + 1][j] + 107.0f/6.0f * zx[i + 2][j] - 13.0f * zx[i + 3][j] + 61.0f/12.0f * zx[i + 4][j] - 5.0f/6.0f * zx[i + 5][j]) / (h*h);
  30.    if (i==n-2)
  31.       return (15.0f/4.0f * zx[i][j] - 77.0f/6.0f * zx[i - 1][j] + 107.0f/6.0f * zx[i - 2][j] - 13.0f * zx[i - 3][j] + 61.0f/12.0f * zx[i - 4][j]) - 5.0f/6.0f * zx[i - 5][j]/ (h*h);
  32.    //if (i == 0 || i == n - 1)
  33.       return 0.0f;
  34. }
  35.  
  36.  
  37. int _tmain(int argc, _TCHAR* argv[])
  38. {
  39. //PARAMETERS
  40.  
  41.    float tmin = 0.0f;
  42.    float tmax = 12.0f;
  43.    int nt = 37000;      //number of points for time
  44.    float dt = (tmax - tmin) / float(nt);
  45.    
  46.    float *time = new float[nt];
  47.    for (int j = 0; j<nt; j++)
  48.       time[j] = tmin + dt*float(j);
  49.  
  50.    float xmin = -30.0f;
  51.    float xmax = 30.0f;
  52.    int nx = 901;        //number of points for x-coordinate
  53.    float *x = new float[nx];
  54.    for (int i = 0; i<nx; i++)
  55.       x[i] = xmin + (xmax - xmin)*float(i) / float(nx);
  56.  
  57.    float v = -0.1f; // velocity of wave
  58.    
  59. //SOLUTION
  60.    float** wave = new float*[nx];       //array of solution
  61.    for (int i = 0; i < nx; i++)
  62.       wave[i] = new float[nt];
  63.  
  64.    float** waved = new float*[nx];      //array of 1st derivative of solution
  65.    for (int i = 0; i < nx; i++)
  66.       waved[i] = new float[nt];
  67.  
  68.    for (int i = 0; i<nx; i++) { //initial conditions
  69.       float xx = x[i];
  70.       wave[i][0] = sin(xx); //wave[i][3] = sin(xx); wave[i][4] = sin(xx);
  71.       waved[i][0] = v*cos(xx); //waved[i][3] = v*cos(xx); waved[i][4] = v*cos(xx);
  72.      
  73.       if (i<6){
  74.         cout<< "i="<< i << "; sinus(" << xx << ")="<< sin(xx) << "; ";
  75.         cout<< "v*cos(xx)=" <<  v*cos(xx) << endl;
  76.       }
  77.  
  78.       for (int j = 1; j<nt; j++){ //initialize rest of memory
  79.         wave[i][j] = 0;
  80.         waved[i][j] = 0;
  81.       }
  82.    }
  83.    cout << endl;
  84.  
  85.    float* wavedd = new float[nx];       //auxiliary variables
  86.    float* waveddot = new float[nx];
  87.  
  88.    for (int i = 0; i<nx; i++) { //j=0;
  89.       wavedd[i] = z2der(nx, i, 0, x, wave);
  90.       waveddot[i] = v*wavedd[i];
  91.     if (i<6){
  92.         cout<< "i="<< i << "; wavedd[" << i <<"] = "<< wavedd[i]  << endl;
  93.       }
  94.    }
  95.    cout << endl;
  96.  
  97.    for (int j = 0; j<nt - 1; j++) {     // Euler method solution
  98.       for (int i = 0; i<nx; i++) {
  99.          waved[i][j + 1] = waved[i][j] + dt * waveddot[i];
  100.          wave[i][j + 1] = wave[i][j] + dt * waved[i][j];
  101.          
  102.          wavedd[i] = z2der(nx, i, j+1, x, wave);
  103.          waveddot[i] = v*wavedd[i];
  104.          
  105.                 if ((i<6)&&(j<6)){
  106.                 cout<< "i="<< i <<", j=" << j << "; wave[" << i <<"," << j <<"] = "<<wave[i][j]<< endl;
  107.         }
  108.       }
  109.    }
  110.    cout << endl;
  111.  
  112.    return 0;
  113. }

 
 
 
 Re: Численное решение волнового уравнения
Сообщение18.08.2015, 22:01 
NT2000 в сообщении #1045526 писал(а):
Думаю ТС уже всё поправил, коль скоро ничего не пишет,
но всё равно опубликую свою версию - одинаково работающую независимо от компилятора.
Добавлено инициация выделенной области, плюс поправки в знаках в граничных точках.
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
  1. // ConsoleApplication1.cpp
  2. //
  3. #include <iostream>
  4. #include <math.h>
  5. #include "stdafx.h"
  6.  
  7. using namespace std;
  8.  
  9. //2nd derivative function
  10. float z2der(int n, int i, int j, float* x, float** zx) {
  11.    float h = x[2] - x[1];
  12.    float zr, zl, zm, zrr, zll;
  13.  
  14.    zm = -30.0f*zx[i][j];
  15.    if (i<n - 1)
  16.       zr = 16.0f*zx[i + 1][j];
  17.    else
  18.       zr = 16.0f;
  19.    if (i>0)
  20.       zl = 16.0f*zx[i - 1][j];
  21.    if (i<n - 2)
  22.       zrr = -zx[i + 2][j];
  23.    if (i>1)
  24.       zll = -zx[i - 2][j];
  25.    
  26.    if (i>1 && i<n-2)
  27.       return (zll + zrr + zl + zr + zm) / (12.0f * h * h);
  28.    if (i == 1)
  29.       return (15.0f/4.0f * zx[i][j] - 77.0f/6.0f * zx[i + 1][j] + 107.0f/6.0f * zx[i + 2][j] - 13.0f * zx[i + 3][j] + 61.0f/12.0f * zx[i + 4][j] - 5.0f/6.0f * zx[i + 5][j]) / (h*h);
  30.    if (i==n-2)
  31.       return (15.0f/4.0f * zx[i][j] - 77.0f/6.0f * zx[i - 1][j] + 107.0f/6.0f * zx[i - 2][j] - 13.0f * zx[i - 3][j] + 61.0f/12.0f * zx[i - 4][j]) - 5.0f/6.0f * zx[i - 5][j]/ (h*h);
  32.    //if (i == 0 || i == n - 1)
  33.       return 0.0f;
  34. }
  35.  
  36.  
  37. int _tmain(int argc, _TCHAR* argv[])
  38. {
  39. //PARAMETERS
  40.  
  41.    float tmin = 0.0f;
  42.    float tmax = 12.0f;
  43.    int nt = 37000;      //number of points for time
  44.    float dt = (tmax - tmin) / float(nt);
  45.    
  46.    float *time = new float[nt];
  47.    for (int j = 0; j<nt; j++)
  48.       time[j] = tmin + dt*float(j);
  49.  
  50.    float xmin = -30.0f;
  51.    float xmax = 30.0f;
  52.    int nx = 901;        //number of points for x-coordinate
  53.    float *x = new float[nx];
  54.    for (int i = 0; i<nx; i++)
  55.       x[i] = xmin + (xmax - xmin)*float(i) / float(nx);
  56.  
  57.    float v = -0.1f; // velocity of wave
  58.    
  59. //SOLUTION
  60.    float** wave = new float*[nx];       //array of solution
  61.    for (int i = 0; i < nx; i++)
  62.       wave[i] = new float[nt];
  63.  
  64.    float** waved = new float*[nx];      //array of 1st derivative of solution
  65.    for (int i = 0; i < nx; i++)
  66.       waved[i] = new float[nt];
  67.  
  68.    for (int i = 0; i<nx; i++) { //initial conditions
  69.       float xx = x[i];
  70.       wave[i][0] = sin(xx); //wave[i][3] = sin(xx); wave[i][4] = sin(xx);
  71.       waved[i][0] = v*cos(xx); //waved[i][3] = v*cos(xx); waved[i][4] = v*cos(xx);
  72.      
  73.       if (i<6){
  74.         cout<< "i="<< i << "; sinus(" << xx << ")="<< sin(xx) << "; ";
  75.         cout<< "v*cos(xx)=" <<  v*cos(xx) << endl;
  76.       }
  77.  
  78.       for (int j = 1; j<nt; j++){ //initialize rest of memory
  79.         wave[i][j] = 0;
  80.         waved[i][j] = 0;
  81.       }
  82.    }
  83.    cout << endl;
  84.  
  85.    float* wavedd = new float[nx];       //auxiliary variables
  86.    float* waveddot = new float[nx];
  87.  
  88.    for (int i = 0; i<nx; i++) { //j=0;
  89.       wavedd[i] = z2der(nx, i, 0, x, wave);
  90.       waveddot[i] = v*wavedd[i];
  91.     if (i<6){
  92.         cout<< "i="<< i << "; wavedd[" << i <<"] = "<< wavedd[i]  << endl;
  93.       }
  94.    }
  95.    cout << endl;
  96.  
  97.    for (int j = 0; j<nt - 1; j++) {     // Euler method solution
  98.       for (int i = 0; i<nx; i++) {
  99.          waved[i][j + 1] = waved[i][j] + dt * waveddot[i];
  100.          wave[i][j + 1] = wave[i][j] + dt * waved[i][j];
  101.          
  102.          wavedd[i] = z2der(nx, i, j+1, x, wave);
  103.          waveddot[i] = v*wavedd[i];
  104.          
  105.                 if ((i<6)&&(j<6)){
  106.                 cout<< "i="<< i <<", j=" << j << "; wave[" << i <<"," << j <<"] = "<<wave[i][j]<< endl;
  107.         }
  108.       }
  109.    }
  110.    cout << endl;
  111.  
  112.    return 0;
  113. }


Да, спасибо. Впредь буду внимательнее. А так действительно, видимо использование "мусора" зависит от компилятора.

 
 
 [ Сообщений: 28 ]  На страницу Пред.  1, 2


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group