implement crank-nicolson in c++

i want to implement crank-nicolson method in c++ according to that :

−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) )

I use gauss-jordan method to solve the system but i can't figure how to implement the above formula.

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;
}

I am assuming that the variable j represents the time steps. In order to implement Crank-Nicolson, you have to pose the problem as a system of linear equations and solve it. The matrix corresponding to the system will be of tridiagonal form, so it is better to use Thomas' algorithm rather than Gauss-Jordan.

The linear system will be of the form A x = b, with x being the vector (..., u(i-1, j+1), u(i, j+1), u(i+1, j+1), ...) and b being the vector (..., r u(i−1,j), 2(1 − r)u(i,j), r u(i+1,j), ...). The i-th row of the matrix A will be of the form (0, ..., 0, −r, 2(1 + r), −r, 0, ..., 0). You will have to be careful with the first and last rows, where you'll have to substitute the boundary conditions for your problem.

A good reference for finite difference methods and Crank-Nicolson in particular is the book by John Strikwerda.

Hope this helps.

C++ is not an equation solving system. = does not have the meaning of equals as it does have in math, but mean assignment.

As such having anything complex like you have on the left side makes no sense, what you probably want to do is to solve the equation so there is a single variable that is being assigned to, possibly with an equation solving program.