Partial Differential Equation of Parabolic Type

KTH ROYAL INSTITUTE OF TECHNOLOGY
School of Engineering Sciences
Course: Applied Numerical Methods

Course code: SF2520
Authors: 1) Andreas Angelou

2) Vasileios Papadimitriou
3) Dionysios Zelios
Lab name: Partial Differential Equation of Parabolic Type
Date of submission: 5/12/2013
CONTENTS
Part A .3
Part B..4
Part C..6
Part D..8
Part E...8
Part F.9
Appendix13
References15
Partial Differential Equation of Parabolic Type

Problem description: We have a metallic rod that is isolated on the right end. At time t=0 we
heat the left end. After some time, the right end will be warmer and then it cools off again.
The density of the rod is given as [kg/m3]. The heat capacity is Cp[J/(kg*C)] and
k[J/(m*s*C)] is the thermal conductivity. The problem can be described by the following:
Length of rod is L[m]. At x=0[m] we have the left end and at x=L[m] we have the
right end.
Initially (t<0) the temperature of the rod is T=0[C].
At time t=0 we heat the left end of the rod with a pulse of temperature T=T0.
The heat pulse lasts for a duration of t=tp[s].
The heat diffusion process can be described by the following partial differential
equation: *Cp*(T/t)=k*(2/x2), t>0, 0<=x<=L
The Boundary Conditions are: 1A) T(0,t)=T0, if 0<=t<=tp
1B) T(0,t)=0, if t>tp
2) /t(L,t)=0
Initial Conditions: 1) T(x,0)=T0, if x=0

2) T(x,0)=0, if 0<x<=L
Part A
Show that with the new values u, and defined by: T=T0*u, x=L* and t=tp*, the problem
can be transformed(scaled) into the following dimensionless form: u/=*(2u/2), >0,
0<<1, with Boundary Conditions:
1a) u(0,)=1, if 0<=<=1
1b) u(0,)=0, if >1
2)u/(1,)=0 and Initial Conditions: 1) u(,0)=1, if =0,
2) u(,0)=0, if 0<<=1.
Show that the only remaining parameter is dimensionless.
Solution
The differential equation of the problem is *Cp*(T/t)=k*(2/x2), t>0, 0<=x<=L.

So we can say that:
t>0 => tp*>0 => >0 (since tp>0)
0<=x<=L => 0<=L*<=L =>(we divide by L) 0<=<=1
(1) =0*u => =0u, (2) x=L* => x=L* and (3) t=tp* => t=tp*
2/x2 = /x(/x )= /(L*)[(0*u)/(L*)] = (0/L2)*( 2u/2).
Thus, the partial differential equation of the problem can be rewritten by using the
relations (1), (2) and (3) as *Cp*(T/t)=k*(2/x2) =>
*Cp*[(T0*u)/(tp*t)]=k*T0/L2*(2u/2) =>(we transfer all constants to the right hand
side) u/t=(k*tp)/(*Cp*L2)*(2u/2). We name =(k*tp)/(*Cp*L2) and we have

u/t=*(2u/2). Also, we can easily see that is dimensionless:
[]=[((J/m*s*C)*s)/((kg/m3)*(J/kg*s)*m2)] => []=[(J*s*m3*kg*C)/(m*s*C*kg*J*m2)] =>
[]=[] (dimensionless).
The Boundary Conditions are: 1A) T(0,t)=T0, if 0<=t<=tp ,
1B) T(0,t)=0, if t>tp
For 0<=t<=tp => 0<=tp*<=tp => 0<=<=1 we have T(0,t)=T0 => T0*u(0,)=0 => u(0,)=1
For t>tp => tp*>tp => >1 we have T(0,t)=0 => T0*u(0,)=0 => u(0,)=0
For x=L we have x=L* => L=L* =>=1 and t=tp*, so /x(L,t)=0 => (T0*u)/(L*)(1,)=0 =>
u/(1,)=0
Hence, the Boundary Conditions can be written as: 1A) u(0,)=1, if 0<=<=1 , 1b) u(0,)=0, if
>1, 2) u/(1,)=0
The Initial Conditions are: 1) T(x,0)=T0, if x=0 ,
2) T(x,0)=0, if 0<x<=L
We have that x=L*, so for x=0 => L*=0 => =0 and T(x,0)=T0 => T0*u(,0)=T0 => u(,0)=1
For x=L => L*=L => =1 and T(x,0)=0 => T0*u(,0)=0 => u(,0)=0.
Hence, the Initial Conditions can be written as 1) u(,0)=1, if =0,
2) u(,0)=0, if 0<<=1
We can consider a=1, for the following parts.
Part B
Discretize the variable using constant stepsize h and central differences to obtain an ODE
system du/dt=A*u+b(), u(0)=u0, where u0 is the zero vector. Show the dimensions and
structures of A, b() and u0. Show also the discretized grid of the -axis you have used and
how the gridpoints are numbered.
Solution
Our model problem is u/=*(2u/2) =>(=1) u/=(2u/2) with Boundary
Conditions u(0,)=1, if 0<=<=1, u(0,)=0, if >1 and u/(1,)=0 and with Initial
Condition u(,0)=1, if =0 and u(,0)=0, if 0<<=1.
We make a discretization with the MoL by introducing an extra gridpoint after the
right boundary point.
Hence i=i*h, where i=0,1,2,,N+1 and h=1/(+1). The ODE system is

ui,k(t,)/=(ui+1,k(,)-2* ui,k(,)+ ui-1,k(,))/ h2, where ui(,) is a time dependent solution
function associated to the space point i. We also know that u(,)= (ui,k+1(,)- ui,k(,))/2 ht
Hence, substituting those partial differential equations to our initial equation and solving
with respect to ui+1,k :
ui+1,k =((2* ht )/ h2 )*(ui+1,k(,)-2* ui,k(,)+ ui-1,k(,)) + ui,k(,)
The Boundary condition can now be transformed:
u(1,)/ =0 => (UN+1- UN-1)/ h2 =0 => UN+1= UN-1
Initial Condition is given u(,0)=1 for =0, and u(,0)=0 for 0<<=1 . The above ODEs can be
written in the matrix form u/=A*u+b() (stiff system of ODEs), u(0)=u0, where u0 is the
zero vector.
2 1
1 2
A=(1/h2)* 0
1

1
1
b()= =(1/h2)*

1

0
0
1
1
0
21
1 2
and
1

0
0
u(0)=
0

0
Part C
We discretize the ODE with the Eulers explicit method using constant stepsize t. In order to
make our calculations efficient we write our code so that the vector Au+b is formed directly.
We store the whole approximate solution (including initial and boundary condition) in a
large matrix U, where the first column is the initial condition while the first and the last row is
the boundary condition.
We use surf to draw a 3D plot of the solution. We make a graph showing a stable solution
and one graph with an unstable solution.
Stable solution
h
0.2
0.004
/h
0.1
Unstable solution
h
0.2
0.013
/h
0.325
Part D
In this part of the lab we compare the explicit method in ode23 and the implicit method in
ode23s, suitable for stiff problems, in the Matlab library. The two functions run under the
similar conditions (same problem, same tolerance) and for three values of the stepsize h
corresponding to N=10, N=20, N=40 grid points on the -axis.
The comparison comprises:
i.
ii.
iii.
The number of time steps needed to reach =2

The cpu-time needed to do each computation
The maximal stepsize that each method could take
The results are presented in the following table:

N
10
20
40
Time steps
Ode23 | Ode23s
458
|
252
1477
|
298
5424
|
347
Cpu-time (sec)
Ode23 | Ode23s
0.562
| 0.807
0.754
| 1.055
1.440
| 1.610
hmax
Ode23
0.0053
0.0006
0.0004
|
|
|
|
Ode23s
0.0456
0.0456
0.0434
We conclude that the number of time steps is considerably smaller when using a stiff
method.
Part E
With the Malab function odeset , options can be set up by the user in order to make the
computation more efficient. By activating the options Jacobian, Jconstant and Jpattern,
we see that we can smooth out the graph for ode-23s and we have the following table:
N
10
20
40
Time steps
Ode23 | Ode23s
458
|
876
1477
|
2099
5424
|
6892
Cpu-time (sec)
Ode23 | Ode23s
0.562
| 0.710
0.754
| 1.019
1.440
| 2.396
hmax
Ode23
0.0053
0.0006
0.0004
|
|
|
|
Ode23s
0.0046
0.0012
0.0003
Part F
We visualize the result of a successful computation from part D in graphs.
The following graphs show two-dimensional plots of u(t,), 0<,=1 at four timepoints:
=0.5, 1, 1.5, 2.
Ode23:
Ode23s
The following graphs show three-dimensional plots of u as a function of and .
Ode23:
Ode23s:
We conclude that if we want efficient numerical calculations for parabolic problems we

should use non-stiff solver (ode-23) for time efficiency. However, after using the Jacobian
matrix properties in the odeset the stiff solver (ode-23s)becomes equally efficient but with
less time efficiency. The structure of the Jacobian is sparse. The type of the linear equation
solver that we used in the implicit ODE method was solver for full matrices.
Appendix
clear all
close all
clc
%=========
===========
Solution of parabolic model problem with the Mol:
%Solution of u_t=u_xx,
MoL.
u(x,0)=0
, u(0,t)=1, u_x(1,t)=0
using the
%================ discretization parameters

=============================
dx=0.2;
%stepsize in x direction
x=0:dx:1;
%inner boundary points
N=1/dx;
%total number of stepsizes
dt=0.013;
% timestep
T=10;
%Total number of timesteps
c=(2*dt)/((dx)^2);
%constant given that a=1
%=============== initialize arrays of boundary and initial conditions

====
u0=[]; b=[];
%initiating initial and boundary
values
u0(1)=1;
%setting the initial value for
=0
for i=2:N+1
%setting the initial values for
>0
u0(i)=0;
end
for j=1:T
%setting the boundary values for
t<1
b0(j)=1;
end
u=u0';
b=b0;
upast = zeros(1,N+1);
U=zeros(N+1,T);
U(:,1)=u;
solution matrix
U(1,:)=b;
solution matrix
% initial value matrix

%upper boundary condition matrix
%initiating the backwards time line
%initiating the solution matrix U
%adding initial condition to the
%adding boundary condition to the
%==================== Creating the rest matrix

==========================
for i=2:T
upast=u;
%condition for saving the backwards
time line
if (dt*i)>1
%setting the boundary values for
t>1
u(1)=0;
end
for j=2:N
%creating the rest of the solution
matrix
u(j)=c*upast(j+1)-(2*c-1)*upast(j)+c*upast(j-1);
U(:,i)=u;
end
u(N+1)=u(N-1);
%setting the lower boundary condition
end
%3-D plot
surf(U)
title('solution of u_t=u_xx, u(x,0)=[1,0,....,0],
u(0,t<1)=[1,1,....,1] u(0,t=1)=0 using the MoL.')
xlabel('space variable x')
ylabel('time variable y')
zlabel('solution variable u(x,t)')
%=========
===========
Solution of parabolic model problem with the Mol:
%Solution of u_t=u_xx, u(x,0)=[1,0,....,0], u(0,t<1)=[1,1,....,1]

u(0,t=1)=0 using the MoL. Discretized system store in heatflow.m
clear all
clc
global N hx
%discretization parameters
N=40;
%number of inner x points
hx=1/(N+1);
%stepsize in x direction
x=hx:hx:1-hx;
%inner points
u0(1)=1;
for i=2:N
u0(i)=0;
% initial values
end
t0=0; tend=2;
%time interval
options = odeset('RelTol',1e-6, 'AbsTol',1e6,'Jacobian',1,'Jconstant',1,'Jpattern',1);
[time,result]=ode23s(@heatflow,[t0,tend],u0,options);
stiff ode-solver
resultp=[zeros(size(time)) result zeros(size(time))];
BCs
xp=0:hx:1;
%add boundary points
mesh(xp,time,resultp)
title('solution of u_t=u_xx, u(x,0)=[1,0,....,0],
u(0,t<1)=[1,1,....,1] u(0,t=1)=0 using the MoL.')
xlabel('space variable x')
ylabel('time variable y')
zlabel('solution variable u(x,t)')
%non%add
%Function used in the parabolic solver for creating the solution

matrix
function rhs=heatflow(t,u)
global N hx
%discretization parameters
rhs=zeros(N,1);
rhs(1)=(-2*u(1)+u(2)+1)/(hx*hx);
rhs(2:N-1)=(u(1:N-2)-2*u(2:N-1)+u(3:N))/(hx*hx);
rhs(N)=(u(N-1)-u(N))/(hx*hx);
end
References
Introduction to computation and modeling for differential equations,
Lennart Edsberg

Partial Differential Equation of Parabolic Type

Uploaded by

Copyright:

Available Formats

Partial Differential Equation of Parabolic Type

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Partial Differential Equation of Parabolic Type

Uploaded by

Copyright:

Available Formats

KTH ROYAL INSTITUTE OF TECHNOLOGY

School of Engineering Sciences

Course: Applied Numerical Methods

Authors: 1) Andreas Angelou

Lab name: Partial Differential Equation of Parabolic Type

Date of submission: 5/12/2013

Partial Differential Equation of Parabolic Type

Initial Conditions: 1) T(x,0)=T0, if x=0

Show that the only remaining parameter is dimensionless.

The differential equation of the problem is *Cp*(T/t)=k*(2/x2), t>0, 0<=x<=L.

side) u/t=(k*tp)/(*Cp*L2)*(2u/2). We name =(k*tp)/(*Cp*L2) and we have

1B) T(0,t)=0, if t>tp

We can consider a=1, for the following parts.

Hence i=i*h, where i=0,1,2,,N+1 and h=1/(+1). The ODE system is

The number of time steps needed to reach =2

The results are presented in the following table:

The following graphs show three-dimensional plots of u as a function of and .

We conclude that if we want efficient numerical calculations for parabolic problems we

Solution of parabolic model problem with the Mol:

%================ discretization parameters

%=============== initialize arrays of boundary and initial conditions

% initial value matrix

%==================== Creating the rest matrix

Solution of parabolic model problem with the Mol:

%Solution of u_t=u_xx, u(x,0)=[1,0,....,0], u(0,t<1)=[1,1,....,1]

%Function used in the parabolic solver for creating the solution

You might also like

The differential equation of the problem is Cp(T/t)=k*(2/x2), t>0, 0<=x<=L.

side) u/t=(ktp)/(CpL2)(2u/2). We name =(ktp)/(Cp*L2) and we have