Re: PDE with RecurrenceTable

• To: mathgroup at smc.vnet.net
• Subject: [mg124832] Re: PDE with RecurrenceTable
• From: Alexei Boulbitch <Alexei.Boulbitch at iee.lu>
• Date: Thu, 9 Feb 2012 05:34:10 -0500 (EST)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com

```Am 03.02.2012 08:13, schrieb Alexei Boulbitch:

> Dear Community,

>

>

>

> I am trying to make a simple numeric FEM solution of parabolic PDE. It seems that RecurrenceTable function is designed exactly for such a job. Indeed, in the Help/RecurrenceTable/Scope/Partial Difference Equations one finds an example. My problem is that this example is only one, and I would say, it is not basic enough.

>

>

>

> My question: do you know some other examples of the use of the RecurrenceTable for this sort of equations??

>

>

>

> I would like to explain: I already went through the MathGroup archive and have seen numerous posts recommending various sophisticated FEM packages. My question is not about them. I want to learn to make simple programs of this sort myself to fast test an equation at hand.

>

>

>

> For example, here is a classical equation of temperature conductivity:

>

>

>

> (\[PartialD]u(x,t))/\[PartialD]t=a^2*(\[PartialD]^2u(x,t))/\[PartialD]x^2

>

>

>

> with a=0.1, the boundary conditions: u[0,t]==1 and u[1,t]==0 and the initial condition u[x,0]== Cos[3 \[Pi]*x/2];

>

> using the explicit method on the rectangular lattice, taken from a textbook:

>

>

>

> u[j,k+1]=\[Sigma]*u[j+1,k]+(1-2*\[Sigma])*u[j,k]+\[Sigma]*u[j-1,k];

>

> \[Sigma]=(a^2*\[Tau])/h^2;

>

>

>

> Tau and h are temporal and spatial step sizes.

>

>

>

> This is the code:

>

>

>

> a = 0.1;

>

> h = 0.1;

>

> \[Tau] = 0.0001;

>

> \[Sigma] = a^2*\[Tau]/h^2;

>

> lst2 = RecurrenceTable[{u[j,

>

>        k + 1] == \[Sigma]*u[j + 1, k] + (1 - 2*\[Sigma])*

>

>         u[j, k] + \[Sigma]*u[j - 1, k], u[j, 0] == Cos[3 Pi*j/20.],

>

>      u[0, k] == 1, u[10, k] == 0}, u, {j, 0, 10}, {k, 0, 100}];

>

>

>

>

>

> Show[{

>

>    ListPlot3D[lst2,

>

>     AxesLabel ->  {Style["t", 16, Italic], Style["x", 16, Italic],

>

>       Style["u", 16, Italic]}, PlotStyle ->  Blue],

>

>    Plot3D[0, {t, 0, 100}, {x, 0, 10},

>

>     PlotStyle ->  {Yellow, Opacity[0.4]}]

>

>    }]

>

>

>

> The solution obtained this way, however, does not show evolution. It is clear that with increasing t the temperature, u, should forget its initial form and approach to a straight line. What is wrong?

>

>

>

> Thank you, Alexei

>

>

>

>

>

> Alexei BOULBITCH, Dr., habil.

>

> IEE S.A.

>

> ZAE Weiergewan,

>

> 11, rue Edmond Reuter,

>

> L-5326 Contern, LUXEMBOURG

>

>

>

> Office phone :  +352-2454-2566

>

> Office fax:       +352-2454-3566

>

> mobile phone:  +49 151 52 40 66 44

>

>

>

> e-mail: alexei.boulbitch at iee.lu<mailto:alexei.boulbitch at iee.lu>

Hi Alexei,

it seems to me that your tau (delta-t) and/or maximum k is far too small

to observe the flattening.

I'll use ListConvolve, because it is much faster; the algorithm is still

the same:

a=1/10; h=1/10; \[Tau]=1/10000;

\[Sigma]=a^2*\[Tau]/h^2;

kmax=121000; (* just to get the zero of u[3/4,t] *)

To avoid memory problems when plotting the result, I'll take ~ 100 data

lines:

lst2 = NestList[

Flatten[{1, ListConvolve[{\[Sigma], 1 - 2*\[Sigma], \[Sigma]}, #1],

0}] & ,

Table[Cos[(3.*Pi*j)/20], {j, 0, 10}],

kmax][[1 ;; All ;; Floor[kmax/100]]];

Now the flattening is visible "to the naked eye":

Show[ListPlot3D[lst2, AxesLabel -> {Style["x", 16, Italic],

Style["t", 16, Italic], Style["u", 16, Italic]},

PlotStyle -> {Lighter[Blue], Specularity[White, 22]},

DataRange -> {{0, 1}, {0, kmax*\[Tau]}}, PlotRange -> All],

Plot3D[0, {x, 0, 1}, {t, 0, kmax*\[Tau]},

PlotStyle -> {Yellow, Opacity[0.4]}]]

Peter

Thank you, Peter,

Works really faster than my code and the boundary conditions are correctly forced.

Thank you once more, Alexei

Alexei BOULBITCH, Dr., habil.

IEE S.A.

ZAE Weiergewan,

11, rue Edmond Reuter,

L-5326 Contern, LUXEMBOURG

Office phone :  +352-2454-2566

Office fax:       +352-2454-3566

mobile phone:  +49 151 52 40 66 44

e-mail: alexei.boulbitch at iee.lu<mailto:alexei.boulbitch at iee.lu>

```

• Prev by Date: Re: recursion won't work for some reason
• Next by Date: I: Re: importing series of file with the same extension
• Previous by thread: Re: PDE with RecurrenceTable
• Next by thread: Find conditions in the coefficients