Re: numerical approximation to the diffusion equation in Mathematica?

• To: mathgroup at smc.vnet.net
• Subject: [mg31957] Re: [mg31945] numerical approximation to the diffusion equation in Mathematica?
• From: Reza Malek-Madani <research at usna.edu>
• Date: Fri, 14 Dec 2001 04:21:20 -0500 (EST)
• Sender: owner-wri-mathgroup at wolfram.com

```I could not figure the equation your code addresses, but if you are
solving nonlinear reaction-diffusion equations, you might be interested in
the following code. It finds the approximate solution to the
initial-boundary value problem (for the ODE version of this problem see
O'Malley's Intro to Singular Pertubation Methods of ODES, p. 170)

u_t = eps u_xx + sin(pi u), u(x,0) = f(x), u(0,t) = -2, u(1,t) = 2

whose steady-state solution has both boundary and internal layers. I use
the method-of-line in t and finite differencing in x. I have a 2D version
of this problem if you are interested.

The last line of the code, which is commented out, saves the animation to
a file.

Reza
************************************
nn = 128; NN = ToString[nn];
eqns = Table[u[i]'[t] == eps/h^2*(u[i + 1][t] - 2 u[i][t] + u[i - 1][t]) +
Sin[Pi*u[i][t]], {i, 1, nn}] /. {u[0][t] -> -2,
u[nn + 1][t] -> 2};
h = 1/(nn + 1);
eps = 0.0001; x = Table[i*h, {i, nn}];
initdata = Table[u[i][0] == 4*x[[i]] - 2, {i, nn}];
eqns = Flatten[Join[eqns, initdata]];
vars = Table[u[i][t], {i, nn}];
sol = NDSolve[eqns, vars, {t, 0, 10}];
graph = Table[ListPlot[Flatten[{-2, Table[First[u[i][t] /. sol], {i, nn}],
2}], PlotJoined -> True, PlotLabel -> StringJoin["N =", NN, ",
eps =", ToString[eps], ", t = ", ToString[t]], Frame -> True,
GridLines -> Automatic, FrameLabel -> "Method of Line",
Background -> GrayLevel[0.5], PlotStyle -> RGBColor[1, 0, 0]],
{t, 0, 10, 0.1}];
(*Export["SingPertFD.gif", graph, "GIF", ImageSize -> {500, 500}]  *)

-------------------------------------------------------------------------
Reza Malek-Madani               Director of Research
Research Office, MS 10m         Phone: 410-293-2504 (FAX -2507)
589 McNair Road                 DSN:      281-2504
U.S. Naval Academy              Nimitz Room 17 in ERC
Annapolis MD 21402-5031         Email: research at usna.edu

--------------------------------------------------------------------------

On Thu, 13 Dec 2001, 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];
> //
>         }
>

```

• Prev by Date: Re: Re: Solve[] for equations?
• Next by Date: Re: Re: Solve[] for equations?
• Previous by thread: numerical approximation to the diffusion equation in Mathematica?
• Next by thread: Re: numerical approximation to the diffusion equation in Mathematica?