"конвертировать" matlab в c ++ - числовые pde -

Question

Я хочу решить pde в соответствии с этим кодом Matlab (часть):

.............
% 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  
............ 

Я сделал это в 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);


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

Проблема в том, что я получаюфайл данных только со значениями 0 1.1 -3.58979e-09.Может быть проблема в индексации? Или я все испортила?

Answer 1

У вас есть несколько серьезных проблем:

1. У вас есть проблемы со всеми для петель.Вы всегда берете одно значение больше, чем существует в массивах.Поскольку вы начинаете свои циклы с нуля, вы должны использовать < вместо <=.

2. Вы используете недопустимые индексы в массивах: u[J][Nt-1] вместо u[i][j] в цикле инициализациии u[J][n] вместо u[J-1][n] в основном цикле.Например, цикл инициализации должен быть следующим:

   for (int i = 0; i < J; i++) {
       for (int j = 0; j < Nt-1; j++) {
           u[i][j] = 0;
       }
   }
   

   Странно, что вы не получаете исключение нарушения доступа, поскольку вы пытаетесь получить доступ к значениям вне массива.

3. Вы проверяете n == 1 вместо n == 0 для первой итерации

4. Чтобы вычисление было равным Matlab, вы должны вычислить t = (n + 1) * dt;, поскольку в Matlab вы n запускаетеот 1.

5. Вместо цикла for(j=0;j<=J-2;j++) должно быть for(j=1; j < J-1; j++)

6. f(x) в первом цикле не имеет смыслаточно так же как gl(xl, t) & gr(xr, t) в начале основного цикла.

Кстати, как вы определили эти функции?Также могут быть проблемы.

Тем не менее, вы должны исправить эти проблемы, прежде чем сможете продолжить тестирование приложения.