"convert" matlab to c++ -- numerical pde's --_问答_开发者

i want to solve a pde according to this matlab code (part of it):

.............
% Evaluate the initial conditions
x = xl : dx : xr;              % generate the grid point
% f(1:J+1) since array index starts from 1
f = sin(pi*x) + sin(2*pi*x);  

% store the solution at all grid points for all time steps
u = zeros(J+1,Nt);   

% Find the approximate solution at each time step
for n = 1:Nt
    t = n*dt;         % current time
    % boundary condition at left side
    gl = sin(pi*xl)*exp(-pi*pi*t)+sin(2*pi*xl)*exp(-4*pi*pi*t);   
    % boundary condition at right side
    gr = sin(pi*xr)*exp(-pi*pi*t)+sin(2*pi*xr)*exp(-4*pi*pi*t);  
    if n==1    % first time step
       for j=2:J    % interior nodes     
       u(j,n) = f(j) + mu*(f(j+1)-2*f(j)+f(j-1));
       end     
       u(1,n) = gl;   % the left-end point
       u(J+1,n) = gr; % the right-end point 
    else 
       for j=2:J    % interior nodes
         u(j,n)=u(j,n-1)+mu*(u(j+1,n-1)-2*u(j,n-1)+u(j-1,n-1));
       end
       u(1,n) = gl;   % the left-end point
       u(J+1,n) = gr; % the right-end point 
    end
 end

% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
surf(x,tt, u');     % 3-D surface plot  
............

I did this in c++ :

double f(double x){

return sin(pi*x) + sin(2.0*pi*x);

}

double gl(double x,double t){

    return sin(pi*x)*exp(-pi*pi*t) +sin(2.0*pi*x)*exp(-4.0*pi*pi*t);
}

double gr(double x,double t){

    return sin(pi*x)*exp(-pi*pi*t) +sin(2.0*pi*x)*exp(-4.0*pi*pi*t);
}

using namespace std;

const int Nt=50;  // number of time steps
const int J=10; // number of division for x
double xl=0,xr=1.0; //left and right boundaries
double tf=0.1; //final simulation time
double dt,dx; //dx is the mesh size
double x;
double u[J][Nt];


int main()
{
dx=(xr-xl)/J;
dt=tf/Nt;
double m=dt/(dx*dx);

开发者_C百科
//Evaluate the initial conditions
//generate the grid point
for (x=xl;x<=xr;x+=dx) {
    f(x);
}

//store the solution at all grid points for all time steps
for (int i=0;i<J;i++){
    for (int j=0;j<Nt-1;j++){
        u[i][j]=0;
    }
}

//find the approximate solution at each time step
int n,t,j;
for (n=0;n<Nt;n++){
    t=(n+1)*dt; //current time
        // boundary condition at left side
        gl(xl,t);
        //boundary condition at right side
        gr(xr,t);

        if (n==0) {//first time step
            for (j=1;j<J-1;j++){
                u[j][n]=f(j)+m*(f(j+1)-2.0*f(j)+f(j-1));
                }
        u[0][n]=gl(xl,t);
        u[J-1][n]=gr(xr,t); }
        else  {
            for(j=1;j<J-1;j++){
                u[j][n]=u[j][n-1]+m*(u[j+1][n-1]-2.0*u[j][n-1]+u[j-1][n-1]);
                }
        u[0][n]=gl(xl,t);
        u[J-1][n]=gr(xr,t);
        }

    }

ofstream data;
data.open("Data.dat");


for (int tt=dt;tt<=Nt*dt;tt+=dt){
    data<< tt <<"\t"<<x<<"\t"<< u[J][Nt]<<endl;

}

cout <<"\nFiles written! "<<endl;

    return 0;
}

The problem is that i get a data file with only values 0 1.1 -3.58979e-09. Maybe the problem is in the indexing?Or i messed up completely?

You have several serious problems here:

You have problems with all for lops. You always takes one value more than exist in arrays. Since you start your loops with zero you should use < instead of <=.
You use invalid indices in arrays: u[J][Nt-1] instead of u[i][j] in the initialization loop, and u[J][n] instead of u[J-1][n] in the main loop. For example the initialization loop should be:
```
for (int i = 0; i < J; i++) {
    for (int j = 0; j < Nt-1; j++) {
        u[i][j] = 0;
    }
}
```
It's strange that you don't get access violation exception, since you are trying to access values outside the array.
You check n == 1 instead of n == 0 for the first iteration
For the calculation to be equal with Matlab you should calculate t = (n + 1) * dt;, since in the Matlab you n starts from 1.
Instead the loop for(j=0;j<=J-2;j++) there should be for(j=1; j < J-1; j++)
The f(x) in the first loop has no sense, just like the gl(xl, t) & gr(xr, t) in the beginning of the main loop.