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