# i want help in mathematica programing

1 view

### Farkhanda

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)

Farkhanda

### jf alcover

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.}}