Re: numerical approximation to the diffusion equation in Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg31971] Re: numerical approximation to the diffusion equation in Mathematica?
- From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
- Date: Fri, 14 Dec 2001 16:52:59 -0500 (EST)
- Organization: Universitaet Leipzig
- References: <9vcdgc$35r$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi,
I realy love this OO style, it's so short !
Since you don't like to write down your
equation and the boundary conditions in a
human readable way here is the version
to solve
D[u[x,t]]==DD*D[u[x,t],{x,2}]
Element[x,Interval[{0,1}]
u[x,0]=If[x<0.5,0,1]
u[0,t]==u[1,t]
So here come the Mathematica input:
(* Dump explicit Euler *)
DiffusionIter[DD_, h_, dx_, u_] :=
u + h*DD/(dx^2)*ListConvolve[{1, -2, 1}, u, 2]
(* Initial concentration *)
u0 = Table[If[i < 11, 0, 1], {i, 20}]
(* Solve the equation with h=0.001, dx=0.0.5 , Diffusion constant =1/2
*)
ll = NestList[DiffusionIter[0.5, 0.001, 0.01, #] &, u0, 50];
(* Draw the result *)
ListPlot3D[ll, Mesh -> False, MeshRange -> {{0, 1}, {0, 0.001*50}}]
Regards
Jens
Bob Marlow wrote:
>
> 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];
> //
> }