реализовать Crank-Николсон в C ++

Question

Я хочу реализовать метод Crank-Nicolson в c ++ согласно этому:

-ru (i − 1, j + 1) + 2 (1 + r) u (i, j +)1) - ru (i + 1, j + 1) = 2 (1 - r) u (i, j) + r (u (i − 1, j) + u (i + 1, j))

Я использую метод Гаусса-Джордана для решения системы, но не могу понять, как реализовать приведенную выше формулу.

const double pi=3.14159265;

double f (double x){
return sin(pi*x);

}

using namespace std;


//gauss-jordan
double* gauss(int n ,double **a){
    double factor;
    double *b,*x;
    x=new double[n];
    b=new double[n];


for (int k=1;k<=n;k++) {
  for (int i=1;i<=n;i++) {
    if (i!=k) {
     factor=a[i][k]/a[k][k];
     for (int j=1; j<=n; j++ ) {
       a[i][j]=a[k][j]*factor-a[i][j];

     }
b[i]=b[k]*factor -b[i];
    }

}
}
for (int i=1;i<=n;i++) {
 x[i]=b[i]/a[i][i];

}
return x;
}

int main()
{
    int i,j,n,m,xd,td;
    double h,k,r;
    double **t,**p;

    //----- initialize double pointer------------------
    double **u;
    u=new double *[m];
    for (int o=1;o<=m;o++){
        u[o]=new double [n];
    }
    //----- end of initialization-----------------------

   cout <<"\nEnter the value of x dimension : \n";
   cin >> xd;
   cout <<"\nEnter the step size h for x dimension : \n";
   cin >>h;
   cout <<"\nEnter the value of time dimension : \n";
   cin >>td;
   cout<<"\nEnter the step k of time dimension : \n";
   cin >>k;
   n=xd/h -1.0;
   m=td/k -1.0;
   cout <<"\n\nThe internal elements of x dimension are :"<<n;
   cout <<"\nThe internal elements of t dimension are :"<<m;
   r=k/(h*h);
   cout <<"\nThe r value is : "<<r;


   //initial conditions
    for (j=0;j<=m;j++){
    u[0][m]=0;
    u[10][m]=0;
    }

    //get the function
    for (i=1;i<n;i++){
    u[i][0]=f(i*h);
    }


    //apply crank-nicolson
    for (i=1;i<n;i++){
        for (j=1;j<n;j++){
              -r*u[i-1][j+1] +2.0*(1.0+r)*u[i][j+1] -r*u[i+1][j+1]=2.0*(1.0-r)*u[i][j] +r*(u[i-1][j]+u[i+1][j]);
    }  // here i can't figure the steps i must follow in order for this to work 



//-----delete double pointer-------------
    for(int o=1;o<m;o++){
    delete [] u[o];
    delete [] u;
    }
//---------------------------------------


    return 0;
}

Answer 1

Я предполагаю, что переменная j представляет временные шаги.Чтобы реализовать Крэнка-Николсона, вы должны представить задачу как систему линейных уравнений и решить ее.Матрица, соответствующая системе, будет иметь трехдиагональную форму, поэтому лучше использовать алгоритм Томаса , а не Гаусса-Джордана.

Линейная система будет иметь вид A x =b, где x - вектор (..., u (i-1, j + 1), u (i, j + 1), u (i + 1, j + 1), ...) и b -вектор (..., ru (i − 1, j), 2 (1 - r) u (i, j), ru (i + 1, j), ...).I-я строка матрицы A будет иметь вид (0, ..., 0, −r, 2 (1 + r), −r, 0, ..., 0).Вам нужно быть осторожным с первой и последней строками, в которых вам придется подставлять граничные условия для вашей задачи.

Хорошей ссылкой для методов конечных разностей и, в частности, Кранка-Николсона, является книга Джона Стрикверда .

Надеюсь, это поможет.

Answer 2

C ++ не является системой решения уравнений. = не имеет значения равно , как в математике, но означает присвоение.

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