Thread: Solving Schrodinger Equation using time-dependent finite method.

I was able to run the basic code using TrueBasic free version. It seems like the plotting commands interface to a runtime library only in that version of BASIC.

To translate the TrueBasic plotting calls to C, you could use a fictional library, "area", which implements the same graph plotting functions as in TrueBasic. A numerical parameter #1 becomes an opaque type AREA*.

tdse.c:

Code:

// PROGRAM tdse
// motion of a wavepacket incident on a potential
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "area.h"

const double PI = 3.14159265358979323846264338327950;
static double Re[1101], Im[1101], Imold[1101];
void parameters(double *x0, double *k0, double *width, double *V0, double *a, 
        double *xmin, double *xmax, double *n, double *dx, double *dx2, 
        double *dt);
void initial_packet(double Re[], double Im[], double x0, double k0,
        double width, double xmin, double n, double dx, double dt);
void set_up_windows(double xmin, double xmax, double V0, 
        AREA **area1, AREA **area2, AREA **area3, AREA **area4);
void draw_potential(double V0, double a, double xmin, double xmax,
        AREA *area4);

int main(void)
{
    void evolve(double Re[], double Im[], double Imold[], double *t, 
            double V0, double a, double dx, double dx2, double dt, 
            double xmin, double n);
    void plots(double Re[], double Im[], double Imold[], double *t, 
            double xmin, double n, double dx, AREA *area1, AREA *area2, 
            AREA *area3, AREA *area4);
    double x0, k0, width, V0, a, xmin, xmax, n, dx, dx2, dt;
    parameters(&x0, &k0, &width, &V0, &a, &xmin, &xmax, &n, &dx, &dx2, &dt);
    initial_packet(Re, Im, x0, k0, width, xmin, n, dx, dt);
    AREA *area1, *area2, *area3, *area4;
    set_up_windows(xmin, xmax, V0, &area1, &area2, &area3, &area4);
    draw_potential(V0, a, xmin, xmax, area4);
    double t = 0.0;
    do {
        evolve(Re, Im, Imold, &t, V0, a, dx, dx2, dt, xmin, n);
        plots(Re, Im, Imold, &t, xmin, n, dx, area1, area2, area3, area4);
    } while (getchar() != EOF);
    closearea(area1);
    closearea(area2);
    closearea(area3);
    closearea(area4);    
    return EXIT_SUCCESS;
}

void parameters(double *x0, double *k0, double *width, double *V0, double *a, 
        double *xmin, double *xmax, double *n, double *dx, double *dx2, 
        double *dt)
{
    *x0 = -15;
    *width = 1;
    *k0 = 2;
    *xmax = 20;
    *xmin = -*xmax;
    *V0 = 2;
    *a = 1;
    *dx = 0.4;
    *dx2 = (*dx)*(*dx);
    *n = (*xmax - *xmin)/(*dx);
    *dt = 0.1;
}

void initial_packet(double Re[], double Im[], double x0, double k0,
        double width, double xmin, double n, double dx, double dt)
    // initial Gaussian wavepacket
{
    double delta2 = width*width;
    double A = pow(2*PI*delta2, -0.25);
    double b = 0.5*k0*dt;
    for (int i=1; i <= n; i++) {
        double x = xmin + (i-1)*dx;
        double arg = 0.25*pow(x - x0, 2)/delta2;
        double e = exp(-arg);
        Re[i] = A*cos(k0*(x - x0))*e;
        Im[i] = A*sin(k0*(x - x0 - 0.5*b))*e;
    }
}

void set_up_windows(double xmin, double xmax, double V0, 
        AREA **area1, AREA **area2, AREA **area3, AREA **area4)
{
    *area1 = openarea(0, 1, 0.75, 1.0);
    setwindow(*area1, xmin, xmax, -0.1, 0.5);
    plotline(*area1, xmin, 0, xmax, 0);
    *area2 = openarea(0, 1, 0.5, 0.75);
    setwindow(*area2, xmin, xmax, -1, 1);
    plotline(*area2, xmin, 0, xmax, 0);
    *area3 = openarea(0, 1, 0.25, 0.5);
    setwindow(*area3, xmin, xmax, -1, 1);
    plotline(*area3, xmin, 0, xmax, 0);
    *area4 = openarea(0, 1, 0, 0.22);
    setwindow(*area4, xmin, xmax, -V0, V0);
}

void draw_potential(double V0, double a, double xmin, double xmax, 
        AREA *area4)
{
    double V(double x, double V0, double a);
    AREA *area = area4;
    setcolor(area, "red");
    plotline(area, xmin, 0, a, 0);
    plotline(area, a, 0, a, V0);
    plotline(area, a, V0, xmax, V0);
    setcolor(area, "black/white");
    printf("total probability = ");
}

void plots(double Re[], double Im[], double Imold[], double *t, double xmin, 
        double n, double dx, AREA *area1, AREA *area2, AREA *area3, AREA *area4)
{
    AREA *area = area2;
    cleararea(area);
    double xmax = 0.0;
    plotline(area, xmin, 0, xmax, 0);
    for (int i=0; i <= n; i++) {
        double x = xmin + (i - 1)*dx;
        plotpoint(area, x, Im[i]);
    }
    area = area1;
    cleararea(area);
    plotline(area, xmin, 0, xmax, 0);
    double Psum = 0;
    for (int i=1 ; i <= n; i++) {
        double x = xmin + (i - 1)*dx;
        double P = Re[i]*Re[i] + Im[i]*Imold[i];
        Psum += P*dx;
        plotpoint(area, x, P);
    }
    area = area4;
    setcursor(area, 1, 20);
    printf("        ");
    setcursor(area, 1, 20);
    printf("%g", Psum);
}

void evolve(double Re[], double Im[], double Imold[], double *t,
        double V0, double a, double dx, double dx2, double dt,
        double xmin, double n)
{
    double V(double x, double V0, double a);
    for (int i=1; i <= n; i++) {
        double x = xmin + (i-1)*dx;
        double HIm = V(x,V0,a)*Im[i] - 0.5*(Im[i+1] - 2*Im[i] + Im[i-1]) 
                / dx2;
        // real part defined at multiples of dt
        Re[i] += HIm*dt;
    }
    for (int i=1; i <= n; i++) {
        double x = xmin + (i-1)*dx;
        // dt/2 earlier than real part
        Imold[i] = Im[i];
        double HRe = V(x,V0,a)*Re[i] - 0.5*(Re[i+1] - 2*Re[i] + Re[i-1])
                / dx2;
        // dt/2 later than real part
        Im[i] -= HRe*dt;
    }
    *t += dt;        // time of real part
}

double V(double x, double V0, double a)
{
    // step potential
    if (x > a) {
        return V0;
    }else{
        return 0;
    }
}

If you want the graphics in C, you'll have to implement a compatible set of calls that manipulate the AREA objects. I don't really know how they're supposed to work, so these are just stub functions to get the thing to compile.

area.h:

Code:

#ifndef AREA_H
#define AREA_H

struct AREA;
typedef struct AREA AREA;

AREA *openarea(double a, double b, double c, double d);
void closearea(AREA *area);
void cleararea(AREA *area);
void setcolor(AREA *area, const char *colorstr);
void setcursor(AREA *area, int row, int col);
void setwindow(AREA *area, double xmin, double xmax, double ymin, 
        double ymax);
void plotline(AREA *area, double x0, double y0, double x1, double y1);
void plotpoint(AREA *area, double x0, double y0);

#endif // AREA_H

area.c:

Code:

#include "area.h"
#include <stdlib.h>
#include <assert.h>

struct AREA {
    double a;
    double b;
    double c;
    double d;
    void *openfn;
};

AREA *openarea(double a, double b, double c, double d)
{
    AREA *area = malloc(sizeof *area);
    assert(area != NULL);
    area->a = a;
    area->b = b;
    area->c = c;
    area->d = d;
    area->openfn = openarea;
    return area;
}

void closearea(AREA *area)
{
    assert(area != NULL);
    assert(area->openfn == openarea);
    free(area);
}

void cleararea(AREA *area) { }
void setcolor(AREA *area, const char *colorstr) { }
void setcursor(AREA *area, int row, int col) { }
void setwindow(AREA *area, double xmin, double xmax, double ymin, double ymax)
{ }
void plotline(AREA *area, double x0, double y0, double x1, double y1) { }
void plotpoint(AREA *area, double x0, double y0) { }

Maybe you replace it with gnuplot calls, but I'm not sure if gnuplot does animation like the TrueBasic runtime calls do.

You shouldn't use assert to test the return value for malloc.

Code:

AREA *openarea(double a, double b, double c, double d)
{
    AREA *area = malloc(sizeof *area);
    assert(area != NULL);
    ...

The reason is that it is only a debuging tool -> When you are finished debugging, you are supposed to define NDEBUG and all of the asserts are pre-processed out of existence.

Originally Posted by c99tutorial

I was able to run the basic code using TrueBasic free version. It seems like the plotting commands interface to a runtime library only in that version of BASIC.

To translate the TrueBasic plotting calls to C, you could use a fictional library, "area", which implements the same graph plotting functions as in TrueBasic. A numerical parameter #1 becomes an opaque type AREA*.

tdse.c:

Code:

// PROGRAM tdse
// motion of a wavepacket incident on a potential
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "area.h"

const double PI = 3.14159265358979323846264338327950;
static double Re[1101], Im[1101], Imold[1101];
void parameters(double *x0, double *k0, double *width, double *V0, double *a, 
        double *xmin, double *xmax, double *n, double *dx, double *dx2, 
        double *dt);
void initial_packet(double Re[], double Im[], double x0, double k0,
        double width, double xmin, double n, double dx, double dt);
void set_up_windows(double xmin, double xmax, double V0, 
        AREA **area1, AREA **area2, AREA **area3, AREA **area4);
void draw_potential(double V0, double a, double xmin, double xmax,
        AREA *area4);

int main(void)
{
    void evolve(double Re[], double Im[], double Imold[], double *t, 
            double V0, double a, double dx, double dx2, double dt, 
            double xmin, double n);
    void plots(double Re[], double Im[], double Imold[], double *t, 
            double xmin, double n, double dx, AREA *area1, AREA *area2, 
            AREA *area3, AREA *area4);
    double x0, k0, width, V0, a, xmin, xmax, n, dx, dx2, dt;
    parameters(&x0, &k0, &width, &V0, &a, &xmin, &xmax, &n, &dx, &dx2, &dt);
    initial_packet(Re, Im, x0, k0, width, xmin, n, dx, dt);
    AREA *area1, *area2, *area3, *area4;
    set_up_windows(xmin, xmax, V0, &area1, &area2, &area3, &area4);
    draw_potential(V0, a, xmin, xmax, area4);
    double t = 0.0;
    do {
        evolve(Re, Im, Imold, &t, V0, a, dx, dx2, dt, xmin, n);
        plots(Re, Im, Imold, &t, xmin, n, dx, area1, area2, area3, area4);
    } while (getchar() != EOF);
    closearea(area1);
    closearea(area2);
    closearea(area3);
    closearea(area4);    
    return EXIT_SUCCESS;
}

void parameters(double *x0, double *k0, double *width, double *V0, double *a, 
        double *xmin, double *xmax, double *n, double *dx, double *dx2, 
        double *dt)
{
    *x0 = -15;
    *width = 1;
    *k0 = 2;
    *xmax = 20;
    *xmin = -*xmax;
    *V0 = 2;
    *a = 1;
    *dx = 0.4;
    *dx2 = (*dx)*(*dx);
    *n = (*xmax - *xmin)/(*dx);
    *dt = 0.1;
}

void initial_packet(double Re[], double Im[], double x0, double k0,
        double width, double xmin, double n, double dx, double dt)
    // initial Gaussian wavepacket
{
    double delta2 = width*width;
    double A = pow(2*PI*delta2, -0.25);
    double b = 0.5*k0*dt;
    for (int i=1; i <= n; i++) {
        double x = xmin + (i-1)*dx;
        double arg = 0.25*pow(x - x0, 2)/delta2;
        double e = exp(-arg);
        Re[i] = A*cos(k0*(x - x0))*e;
        Im[i] = A*sin(k0*(x - x0 - 0.5*b))*e;
    }
}

void set_up_windows(double xmin, double xmax, double V0, 
        AREA **area1, AREA **area2, AREA **area3, AREA **area4)
{
    *area1 = openarea(0, 1, 0.75, 1.0);
    setwindow(*area1, xmin, xmax, -0.1, 0.5);
    plotline(*area1, xmin, 0, xmax, 0);
    *area2 = openarea(0, 1, 0.5, 0.75);
    setwindow(*area2, xmin, xmax, -1, 1);
    plotline(*area2, xmin, 0, xmax, 0);
    *area3 = openarea(0, 1, 0.25, 0.5);
    setwindow(*area3, xmin, xmax, -1, 1);
    plotline(*area3, xmin, 0, xmax, 0);
    *area4 = openarea(0, 1, 0, 0.22);
    setwindow(*area4, xmin, xmax, -V0, V0);
}

void draw_potential(double V0, double a, double xmin, double xmax, 
        AREA *area4)
{
    double V(double x, double V0, double a);
    AREA *area = area4;
    setcolor(area, "red");
    plotline(area, xmin, 0, a, 0);
    plotline(area, a, 0, a, V0);
    plotline(area, a, V0, xmax, V0);
    setcolor(area, "black/white");
    printf("total probability = ");
}

void plots(double Re[], double Im[], double Imold[], double *t, double xmin, 
        double n, double dx, AREA *area1, AREA *area2, AREA *area3, AREA *area4)
{
    AREA *area = area2;
    cleararea(area);
    double xmax = 0.0;
    plotline(area, xmin, 0, xmax, 0);
    for (int i=0; i <= n; i++) {
        double x = xmin + (i - 1)*dx;
        plotpoint(area, x, Im[i]);
    }
    area = area1;
    cleararea(area);
    plotline(area, xmin, 0, xmax, 0);
    double Psum = 0;
    for (int i=1 ; i <= n; i++) {
        double x = xmin + (i - 1)*dx;
        double P = Re[i]*Re[i] + Im[i]*Imold[i];
        Psum += P*dx;
        plotpoint(area, x, P);
    }
    area = area4;
    setcursor(area, 1, 20);
    printf("        ");
    setcursor(area, 1, 20);
    printf("%g", Psum);
}

void evolve(double Re[], double Im[], double Imold[], double *t,
        double V0, double a, double dx, double dx2, double dt,
        double xmin, double n)
{
    double V(double x, double V0, double a);
    for (int i=1; i <= n; i++) {
        double x = xmin + (i-1)*dx;
        double HIm = V(x,V0,a)*Im[i] - 0.5*(Im[i+1] - 2*Im[i] + Im[i-1]) 
                / dx2;
        // real part defined at multiples of dt
        Re[i] += HIm*dt;
    }
    for (int i=1; i <= n; i++) {
        double x = xmin + (i-1)*dx;
        // dt/2 earlier than real part
        Imold[i] = Im[i];
        double HRe = V(x,V0,a)*Re[i] - 0.5*(Re[i+1] - 2*Re[i] + Re[i-1])
                / dx2;
        // dt/2 later than real part
        Im[i] -= HRe*dt;
    }
    *t += dt;        // time of real part
}

double V(double x, double V0, double a)
{
    // step potential
    if (x > a) {
        return V0;
    }else{
        return 0;
    }
}

If you want the graphics in C, you'll have to implement a compatible set of calls that manipulate the AREA objects. I don't really know how they're supposed to work, so these are just stub functions to get the thing to compile.

area.h:

Code:

#ifndef AREA_H
#define AREA_H

struct AREA;
typedef struct AREA AREA;

AREA *openarea(double a, double b, double c, double d);
void closearea(AREA *area);
void cleararea(AREA *area);
void setcolor(AREA *area, const char *colorstr);
void setcursor(AREA *area, int row, int col);
void setwindow(AREA *area, double xmin, double xmax, double ymin, 
        double ymax);
void plotline(AREA *area, double x0, double y0, double x1, double y1);
void plotpoint(AREA *area, double x0, double y0);

#endif // AREA_H

area.c:

Code:

#include "area.h"
#include <stdlib.h>
#include <assert.h>

struct AREA {
    double a;
    double b;
    double c;
    double d;
    void *openfn;
};

AREA *openarea(double a, double b, double c, double d)
{
    AREA *area = malloc(sizeof *area);
    assert(area != NULL);
    area->a = a;
    area->b = b;
    area->c = c;
    area->d = d;
    area->openfn = openarea;
    return area;
}

void closearea(AREA *area)
{
    assert(area != NULL);
    assert(area->openfn == openarea);
    free(area);
}

void cleararea(AREA *area) { }
void setcolor(AREA *area, const char *colorstr) { }
void setcursor(AREA *area, int row, int col) { }
void setwindow(AREA *area, double xmin, double xmax, double ymin, double ymax)
{ }
void plotline(AREA *area, double x0, double y0, double x1, double y1) { }
void plotpoint(AREA *area, double x0, double y0) { }

Maybe you replace it with gnuplot calls, but I'm not sure if gnuplot does animation like the TrueBasic runtime calls do.

Whoa! Calm down! I'm just a newbie, as you can see from my little random number generator program. Ok, let's see if I understand the whole thing: you translated the entire program to C; then, you took the functions and library calls specifics to BASIC and replaced them with a dummy library AREA. I really appreciate your work here, but this is what I asked about in the first place, quoting myself: "The only code (that solves it) that I found so far is in BASIC, and I can't use it as a template because it uses some BASIC-only libraries, and I don't know their counterparts in C."

Nonetheless, your code is so much cleaner than what I wrote...

Thank you, not because you solved my former question, but because you were very prestative (don't know if this word exists in english (not sure if there is this word - exists - in english)).

Then again, some of you guys know which library should I use to replace those BASIC specific calls?

Also, it would do if I could replace those library calls with printfs to a file, instead of plotting it directly, it would print the resulting numbers in a file, and I could tell gnuplot to use them... I guess...

Originally Posted by terminusest

Ok, let's see if I understand the whole thing: you translated the entire program to C; then, you took the functions and library calls specifics to BASIC and replaced them with a dummy library AREA.

...

Nonetheless, your code is so much cleaner than what I wrote...

It was my intention to make it exactly the same just to bridge the language gap. Didn't you write the BASIC code? Whoever wrote it must know how the method works -- and that somebody is definitely not me. In any case the C version gave me the same numerical answer as the BASIC version.

Originally Posted by terminusest

Also, it would do if I could replace those library calls with printfs to a file, instead of plotting it directly, it would print the resulting numbers in a file, and I could tell gnuplot to use them... I guess...

Yes, then just fill in the AREA functions with printf's that gnuplot will understand it. The only difference I can see is that the BASIC version is using a "clear" function to animate the graph, and I think gnuplot is intended for non-moving plots, but I could be wrong. For example WINDOW #2 followed by CLEAR must involve only the "#2" area, so the C translation becomes

Code:

 area = area2; // WINDOW #2
 cleararea(area);  // CLEAR

If you want the graphing done the same way, best would be to sit down with gnuplot and TrueBASIC side by side and look for the appropriate "translations" for each plotting command to accomplish the same behaviour. For example how do you partition of an "area" of the window in gnuplot, how do you "clear" it, how do you draw a line from x0,y0 to x1,y1, etc... if you know these things, then the AREA functions could in principle be used to emulate this functionality of TrueBASIC.

The only difference I can see is that the BASIC version is using a "clear" function to animate the graph, and I think gnuplot is intended for non-moving plots, but I could be wrong.

Well, in fact it does not "animate" in the sense that it refreshes the graph 30 times per second, but as you can see on my little number generator program some posts above that I can clear the graph every half second and plot a new graph.

If you want the graphing done the same way, best would be to sit down with gnuplot and TrueBASIC side by side and look for the appropriate "translations" for each plotting command to accomplish the same behaviour. For example how do you partition of an "area" of the window in gnuplot, how do you "clear" it, how do you draw a line from x0,y0 to x1,y1, etc... if you know these things, then the AREA functions could in principle be used to emulate this functionality of TrueBASIC.

I think my limited understanding of programming is what is hindering me from translating the whole thing. I do not fully understand how those functions with multiple arguments - like void parameters(double *x0, double *k0, double *width, double *V0, double *a, double *xmin, double *xmax, double *n, double *dx, double *dx2, double *dt); - really work. I'm studying the code right now, and I'm trying to grasp how can I open a screen on X11 with just C libraries.

Any clues?

EDIT: Ok, after some googling I found it's impossible to really "draw" something on console, unless you are after some ASCII art... What about OpenGL, then?

Advices?

Your first question seems like its about pointer notation in C. Notice that in the BASIC version you have a function like this

SUB parameters(x0,k0,width,V0,a,xmin,xmax,n,dx,dx2,dt)
...
END SUB

Notice that inside the function, the code is allowed to modify the variables. To accomplish the equivalent in C one must use pointers. Basically in practice it just means you prefix everything with * in the parameters list

void parameters(double *x0, ...);

And then you must prefix the passed-in local variables with & when you actually call the function.

parameters(&x0, ...);

When dealing with pointers many references tend to go into too complicated explanations and start talking about unnecessary things like function pointers and dynamic memory allocation, linked lists, etc., before the basics are covered. I recommend reading this online chapter and trying the examples out to get the hang of it: Practical C Programming, 3rd Edition -- Sample chapter

Your second question is about drawing primitives. If you're on X11 I suppose you could directly use the X api, but I don't have experience with this and from what I hear it is very painful to use. Gtk has a much better API and is pretty much standard on Linux, and if you include the DLL in your program should work on Windows as well. here is the chapter on drawing primitives: Drawing Primitives

If you want to use this I would start with a simple example first like make a program which draws a simple plot or something to get the hang of it.