[Date Index] [Thread Index] [Author Index]

numerical approximation to the diffusion equation in Mathematica?

• To: mathgroup at smc.vnet.net
• Subject: [mg31945] numerical approximation to the diffusion equation in Mathematica?
• From: bobmarlow at postmaster.co.uk (Bob Marlow)
• Date: Thu, 13 Dec 2001 01:08:53 -0500 (EST)
• Sender: owner-wri-mathgroup at wolfram.com

I have written a very simple numerical
approximation to the diffusion equation -
C++ source included below.

Do you know how I would go about repeating
this in Mathematica, especially drawing a graph
of the results? The difference equation
is from Strauss.

(Excuse all the code comments, I'm doing
some multi-tasking here, i.e. as well as
learning about Inverse Problems (Basic 
PDE study), I'm also learning the OO
aspect of C++. See http://www.inverse-problems.com
for more on that!)


Thanks.

Bob

        // ood //
//
// Macros
//
        #include <conio.h>              // Need for 'getch'
        #include <stdio.h>              // Best for scanf etc
        #include <string.h>             // Best for strcpy etc
        #include <iostream.h>           // cin,cout
        #include <stdlib.h>             // For atoi
        #include "ip.h"                 // Local User-defined Macros
        #include "oodLayer.h"           // Local User-defined Macros

        main(int argc,char *argv[3])
        {
        
//***************************************************************
//                                                              *
//      Polymorphism. Various types. 2 already done:            *
//                                                              *
//      (i)  If you have i=x, you are implicitly forcing x to   *
//           be integer type, or rather, you're relying on the  *
//           compiler to do that.                               *
//                                                              *
//      (ii) By putting i= (int) x, you are explicitly          *
//           converting x to integer type. May be referred to   *
//           as casting.                                        *
//                                                              *
//      We look at another type here:-                          *
//                                                              *
//           For (object) variables, allowing their class to be *
//           determined at run-time. This is called dynamic, or *
//           late binding, or run-time polymorphism.            *
//                                                              *
//                                              25.11.2001.     *
//                                                              *
//***************************************************************

        int j,m,n,obj_type=0,toggle=0;
        Layer* Layers[2];           //Generic layer class pointer
        double* u;
        double s,phi15[]={1,4,6,9,11,14,17,20,17,14,11,9,6,4,1};
        double phi40[]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,\
        20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1};
        s=atof(argv[1]);
        n=atoi(argv[2]);
//
//      Is our layer the oob/ooc 15 points, or the ood 40?
//      We'll let 3rd parameter decide, if set:-
//
        if (argc > 3) obj_type=atoi(argv[3]);
        if (obj_type==40) {
           Layers[0]=new Layer40;
           Layers[1]=new Layer40;
           Layers[0]->array_init(phi40);
           m=40;
        }
        else {
           Layers[0]=new Layer15;
           Layers[1]=new Layer15;
           Layers[0]->array_init(phi15);
           m=15;
        }
//
//      Report 1st layer values:-
//
        u=Layers[0]->fetch();
        cout << "Finite Differences for Diffusion equation" << endl;
        cout << "Initial Data is: " << endl;
        for (j=0; j<m; j++) printf("%4.1f ",*(u+j));
//
//      Now set 2nd up from the first:-
//
        Layers[1]->array_set(s,Layers[0]->fetch());
        u=Layers[1]->fetch();
        cout << endl;
        for (j=0; j<m; j++) printf("%4.1f ",*(u+j));
//
//      Now do n more iterations:-
//
        for (int k=0; k < n; k++) {
           Layers[toggle]->array_set(s,Layers[1-toggle]->fetch());
           u=Layers[toggle]->fetch();
           cout << endl;
           for (j=0; j<m; j++) printf("%4.1f ",*(u+j));
           toggle=1-toggle;
        } //End loop
        }


//
//      Layer class for ood: Polymorphism (dynamic/run-time/late binding).
//
//      Generic layer class:-
//
        class Layer {
protected:
        double points[100];      //Set to maximum here
public:  
//
//      Headers only (MUST have "virtual"), each sub-class has def's :-
//
        virtual double* fetch();
        virtual void array_init(double x[]);
        virtual void array_set(double s,double x[]);
        };
//
//      Need to put dummy bodies in to stop linker
//      errors. Warnings may still appear, but these
//      won't matter, as the member functions will work:-
//
        double* Layer::fetch() {}
        void Layer::array_init(double x[]) {}
        void Layer::array_set(double s,double x[]) {}
//
//      15-point layer:-
//
        class Layer15 : public Layer {
//
public:  
//
//      Method headers:-
//
        double* fetch();
        void array_init(double x[]);
        void array_set(double s,double x[]);
        };
//
//      Object methods:-
//
        double* Layer15::fetch()
        {
                return &points[0];
        }
//
        void Layer15::array_init(double x[])
        {
        int j;
        for (j=0; j<15 ;j++) points[j]=x[j];
        }
//
        void Layer15::array_set(double s,double x[])
        {
        int j;
//
//      For start and end of array, we assume zero:-
//
        points[0]=s*x[1]+(1-2*s)*x[0];
        points[14]=s*x[13]+(1-2*s)*x[14];
//
//      Standard formula for the rest:-
//
        for (j=1; j<14 ;j++) points[j]=s*(x[j+1]+x[j-1])+(1-2*s)*x[j];
//
        }
//
//      Bigger layer:-
//
        class Layer40 : public Layer {
public:  
//
//      Method headers:-
//
        double* fetch();
        void array_init(double x[]);
        void array_set(double s,double x[]);
        };
//
//      Object methods:-
//
        double* Layer40::fetch()
        {
                return &points[0];
        }
//
        void Layer40::array_init(double x[])
        {
        int j;
        for (j=0; j<40 ;j++) points[j]=x[j];
        }
//
        void Layer40::array_set(double s,double x[])
        {
        int j;
//
//      For start and end of array, we assume zero:-
//
        points[0]=s*x[1]+(1-2*s)*x[0];
        points[39]=s*x[38]+(1-2*s)*x[39];
//
//      Standard formula for the rest:-
//
        for (j=1; j<39 ;j++) points[j]=s*(x[j+1]+x[j-1])+(1-2*s)*x[j];
//
        }