1 view

Skip to first unread message

Jul 25, 2003, 11:18:49 AM7/25/03

to

solve the heat equation using explicit difference approximation

with the help of programming in mathematica

let u= u(x,t)

with the given boundary conditions u[0, t] =u[1, t] = 0 for 0 <= t <= 1

and the initial condition u(x,0)=x when 0<=x<=1/2

u(x,0)=1-x when 1/2 <=x<=1

The explicit difference scheme is given as

[u(x,t+delta t)-u(x,t)]/(delta t) =[u(x+delta x,t)-2u(x,t)+u(x-delta x,t)]/(delta x )^2

with a condition

(delta t) / (delta x )^2 <=(1/2)

your sincerely

Farkhanda

Aug 1, 2003, 1:37:13 AM8/1/03

to

Farkhanda <farkhan...@yahoo.com> wrote in message news:<bfrhop$qsr$1...@smc.vnet.net>...

Here is an example (not thoroughly tested,

not the best way, just my way ! )

of what could be done with mma :

The subscripts used are k for x[k] and n for t[n].

ClearAll[u, x, t];

deltat = 0.01;

deltax = 0.2;

(* we check the stability condition : *)

(deltat) / (deltax )^2

0.25

kmax = Floor[1/deltax]

5

nmax = Floor[1/deltat]

100

(* boundary conditions : *)

u[0, n_Integer /; 0 <= n <= nmax] = 0;

u[kmax, n_Integer /; 0 <= n <= nmax] = 0;

u[k_Integer /; 0 <= k <= kmax/2, 0] := k deltax;

u[k_Integer /; kmax/2 <= k <= kmax, 0] := 1 - k deltax;

(* difference scheme : *)

eq = (u[k, n + 1] - u[k, n])/(deltat) == (u[k + 1, n] - 2u[k, n] +

u[k - 1, n])/(deltax )^2;

(* solve for u[k,n+1] : *)

sol = Solve[eq, u[k, n + 1]]

{{u[k, 1 + n] ->

0.01 (100. u[k, n] + 25. (u[-1 + k, n] - 2. u[k, n] + u[1 + k, n]))}}

(* we need another subscript because on left hand side n+1 is not allowed : *)

sol2 = sol[[1]] /. n -> m - 1

{u[k, m] ->

0.01 (100. u[k, -1 + m] +

25. (u[-1 + k, -1 + m] - 2. u[k, -1 + m] + u[1 + k, -1 + m]))}

sol2[[1, 2]] // Simplify

0.25 u[-1 + k, -1 + m] + 0.5 u[k, -1 + m] + 0.25 u[1 + k, -1 + m]

(* we plug this formula into u[k,m] : *)

u[k_Integer, m_Integer /; m < 0] = 0;

u[k_Integer /; k < 0 || k > kmax] = 0;

u[k_Integer /; -kmax <= k <= 2kmax, m_Integer /; m > 0] :=

(* cache is useful for speed *)

u[k, m] = 0.25 u[-1 + k, -1 + m] + 0.50 u[k, -1 + m] + 0.25 u[1 + k, -1 + m];

(* results for the 5 first layers : *)

Table[u[k, n], {k, 0, kmax}, {n, 0, 5}]

{{0., 0., 0., 0., 0., 0.},

{0.2, 0.2, 0.1875, 0.171875, 0.15625, 0.141602},

{0.4, 0.35, 0.3125, 0.28125, 0.253906, 0.229492},

{0.4, 0.35, 0.3125, 0.28125, 0.253906, 0.229492},

{0.2, 0.2, 0.1875, 0.171875, 0.15625, 0.141602},

{0., 0., 0., 0., 0., 0.}}

I hope this will help you to build your own program.

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu