Блог

У меня проблема с программой, которую я пишу на C, которая решает 1D уравнение линейной конвекции. По сути, я инициализировал два массива. Первый массив (u0_array) представляет собой массив единиц с элементами arrays, равными двум за интервал или 0,5 < x < 1. Второй массив (usol_array) служит временным массивом, в котором будет сохранен результат или решение.

Проблема, с которой я сталкиваюсь, — это вложенный цикл for в конце моего кода. Этот цикл применяет уравнение обновления, необходимое для вычисления следующей точки, к каждому элементу массива. Когда я запускаю свой скрипт и пытаюсь распечатать результаты, я получаю всего -1.IND00 для каждой итерации цикла. (Я следую псевдокоду в стиле Matlab, который я также прикрепил ниже) Я очень новичок в C, поэтому моя неопытность показывает. Я не знаю, почему это происходит. Если кто-нибудь может предложить возможное исправление этого, я был бы очень благодарен. Я прикрепил свой код ниже с несколькими комментариями, чтобы вы могли следить за моим мыслительным процессом. Я также прикрепил псевдокод в стиле Matlab, за которым я следую.

 # include <math.h>
# include <stdlib.h>
# include <stdio.h>
# include <time.h>


int main ( )
{
    //eqyation to be solved 1D linear convection -- du/dt   c(du/dx) = 0
    //initial conditions u(x,0) = u0(x)

    //after discretisation using forward difference time and backward differemce space
    //update equation becomes u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]);


    int nx = 41; //num of grid points

    int dx = 2 / (nx - 1); //magnitude of the spacing between grid points

    int nt = 25;//nt is the number of timesteps

    double dt = 0.25; //the amount of time each timestep covers

    int c = 1;  //assume wavespeed


    //set up our initial conditions. The initial velocity u_0
    //is 2 across the interval of 0.5 <x < 1 and u_0 = 1 everywhere else.
    //we will define an array of ones
    double* u0_array = (double*)calloc(nx, sizeof(double));

    for (int i = 0; i < nx; i  )
    {
        u0_array[i] = 1;
    }

    // set u = 2 between 0.5 and 1 as per initial conditions
    //note 0.5/dx = 10, 1/dx 1 = 21
    for (int i = 10; i < 21; i  )
    {
        u0_array[i] = 2;
        //printf("%f, ", u0_array[i]);
    }

    //make a temporary array that allows us to store the solution
    double* usol_array = (double*)calloc(nx, sizeof(double));


    //apply numerical scheme forward difference in
    //time an backward difference in space
    for (int i = 0; i < nt; i  )
    {
        //first loop iterates through each time step
        usol_array[i] = u0_array[i];
        //printf("%f", usol_array[i]);
        

        //MY CODE WORKS FINE AS I WANT UP TO THIS LOOP
        //second array iterates through each grid point for that time step and applies
        //the update equation
        for (int j = 1; j < nx - 1; j  )
        {
            u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
            printf("%f, ", u0_array[j]);
        }
    }

    return EXIT_SUCCESS;
}
 

Для справки также псевдокод, за которым я следую, прикреплен ниже
псевдокода линейной конвекции 1D (стиль Matlab)

Вместо целочисленного деления используйте FP math.
Это позволяет избежать последующего деления на 0 и -1.#IND00

 // int dx = 2 / (nx - 1);   quotient is 0.
double dx = 2.0 / (nx - 1);
 

Код OP не соответствует комментарию

 // u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]
u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
 

Легче увидеть _array , удалена ли избыточность.

 //                                     v correct?
// u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]
u0[j] = usol[j] - c * dt/dx * (usol[j] - usol[j - 1]);
//                                       ^^^^ Correct?
 

То, что вы, вероятно, хотите, это в терминах matlab

 for i = 1:nt
    usol = u0
    u0(2:nx) = usol(2:nx) - c*dt/dx*(usol(2:nx)-usol(1:nx-1))
end%for
 

Это означает, что у вас есть 2 внутренних цикла над пространственным измерением, по одному для каждой из двух векторных операций. Эти два отдельных цикла должны быть явными в коде C

     //apply numerical scheme forward difference in
    //time an backward difference in space
    for (int i = 0; i < nt; i  ) {
        //first loop iterates through each time step
        for (int j = 0; j < nx; j  ) {
            usol_array[j] = u0_array[j];
            //printf("%f", usol_array[i]);
        } 
        

        //second loop iterates through each grid point for that time step 
        //and applies the update equation
        for (int j = 1; j < nx; j  )
        {
            u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
            printf("%f, ", u0_array[j]);
        }
    }