guidance in mathematica programming in pde
- To: mathgroup at smc.vnet.net
- Subject: [mg43628] guidance in mathematica programming in pde
- From: farkhanda_yusaf at hotmail.com (farkhanda)
- Date: Fri, 26 Sep 2003 04:45:40 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
respected sir i have just started the mathematica and facing lot of
problems. during my search on internet I found your group email
address. you repliled to one student like me about the problem in
mathematica.I am sending the problem with my approach of
programming,if possible point out me my mistake and at the same time
if write the programm of this problem again,it will better for me in
the sense that it give me new idea to think on the same line.
waiting for a very positive response.thank you very much
PROBLEM "solve the heat equation using explicit difference
approximation with the help of programming in mathematica"
let u= u(x,t) and heat equation is
derivative of u(x,t) wrt t =double derivative of u(x,t) wrt x
with the given boundary conditions
u(0, t) =0, u(1, t) = 0 for t greater than and equal to 0
and the initial condition u(x,0)=x when x belongs to [0,1/2]
u(x,0)=1-x when x belongs to [1/2,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 is les than equal to (1/2)
my requirement is I need the graph in 3 dimension of u, x and t.
MY APPROACH OF PROGRAMING
ClearAll[u, x, t];
deltat = 0.01;
deltax = 0.2;
(*we check the stability condition :*)
Print["stability test = ", (deltat)/(deltax)^2 <= 1/2];
kmax = Floor[1/deltax];
nmax = Floor[1/deltat];
(*boundary conditions :*)
u[k_Integer /; k <= 0 || k >= kmax, n_Integer] = 0;
(*initial conditions :*)
u[k_Integer /; 0 <= k <= kmax/2, 0] := k deltax;
u[k_Integer /; kmax/2 <= k <= kmax, 0] := 1 - k deltax;
(*difference scheme :*)
(u[k, n + 1] - u[k, n])/(deltat) == (u[k + 1, n] - 2u[k, n] +
u[k - 1, n])/(deltax)^2;
eq = 0.25 u[-1 + k, n] + 0.50 u[k, n] + 0.25 u[1 + k, n]
u[k_Integer /; 0 <= k <= kmax, n_Integer /; n > 0] :=
u[k, n] =
0.25 u[-1 + k, n - 1] + 0.50 u[k, n - 1] + 0.25 u[1 + k, n - 1];
abc = Table[u[k, n], {k, 0, kmax}, {n, 0, nmax - 1}] // Transpose //
MatrixForm
ListPlot3D[abc]