A..10 9_Diffusion_1D.c

9_Diffusion_1D.c

#include<stdio.h>
#include<glsc3d_3.h>

#define N (200)

int main()
{
    double  x;
    double  L = 5;
    double  dx = L / N;
    double  t = 0.0;
    double  dt = 0.0001;
    double  u[N + 2];
    double  tmpu[N + 2];
    double  D = 1.0;
    int     INTV = 100;
    g_init("Window", 600, 320); //Pixel Size
    g_def_scale_2D(0,                   //ID
                   -L*0.5, L*0.5,       //xmin,xmax
                   -0.2, 1.2,     //ymin,ymax
                   20.0, 20.0,          //Window (Left, Top) Position
                   560, 280);           //Window Size (x,y)
    //Set initial calculation
    {
        //Step1
        for(int i = 1;i < N + 1; i ++)
        {
            x = i*dx - L*0.5;
            u[i] = exp(-(x*x));
        }
        //Step2
        u[0] = u[1];
        u[N + 1] = u[N];
    }
    //////////// Start time loop ////////////
    for (int i_time = 0;t < 10 ; i_time++) {
        
        t = i_time * dt;
        
        //////////// Draw Part ////////////
        if (i_time%INTV == 0) {
            g_cls();                    //Clear window
            g_sel_scale(0);             //Select Virtual scale
            g_line_color(0.0, 0.0, 0.0, 1.0);
            g_boundary();               //Draw Boundary
            g_line_width(1.0);
            g_line_color(0.0, 0.0, 0.0, 0.5);
            g_move_2D(-L*0.5, 0.0); g_plot_2D(L*0.5, 0.0);  //X Axis
            g_move_2D(0.0, -1.2); g_plot_2D(0.0, 1.2);      //Y Axis
            //Profile of u
            g_line_width(2.0);
            g_line_color(1.0, 0.0, 0.0, 1.0);
            for (int i=0; i<N-1; i++) {
                g_move_2D(i*dx-L*0.5, u[i]);
                g_plot_2D((i+1)*dx-L*0.5, u[i+1]);
            }
            
            g_text_size(18);
            g_text_color(1.0, 0.0, 0.0, 1.0);
            g_text_2D_virtual(L*0.5 - 0.4, 1.0, "u: -");
            
            g_finish();                 //flush Draw buffer
            g_sleep(0.0);               //Sleep 0.0 sec
        }
        //////////// Calculation Part ////////////
        {
            //Step3
            for(int i = 1;i < N + 1; i ++)
            {
                tmpu[i] = u[i] + D * (u[i - 1] -2 * u[i] + u[i + 1]) / (dx * dx) * dt;
            }
            for(int i = 1;i < N + 1; i ++)
            {
                u[i] = tmpu[i];
            }
            //Step4
            u[0] = u[1];
            u[N + 1] = u[N];
        }
    }
    return 0;
}

• 空間1次元の拡散方程式のグラフを描くプログラム。