Transcript 11.MethodOfLines - Electrical and Computer Engineering

MATH 212 Advanced Calculus 2 for Electrical Engineering
Advanced Calculus 2 for Nanotechnology Engineering
NE 217
Method of Lines
Douglas Wilhelm Harder
Department of Electrical and Computer Engineering
University of Waterloo
Waterloo, Ontario, Canada
Copyright © 2011 by Douglas Wilhelm Harder. All rights reserved.
Method of Lines
Outline
This topic introduces the method of lines
–
–
–
–
Used for solving heat-conduction/diffusion and wave equations
Finite differences discretizes both time and space
Discretize only space
Use a high-quality ODE solver to find the solution over time
2
Method of Lines
Outcomes Based Learning Objectives
By the end of this laboratory, you will:
– You will understand the method of lines
– You will be able to implement it using the Matlab functions implemented
in previous terms
3
Method of Lines
Integration-by-Parts in Higher Dimensions
Up to this point, we have discretized both the space and time
dimensions
– We will look at another approach that discretizes only the space
dimensions
4
Method of Lines
Discretizing in Space
Consider the heat-conduction/diffusion equation
1 u
  2u
 t
Discretizing the space component in one dimension gives:


u  x, t   2  u  x  h, t   2u  x, t   u  x  h, t  
t
h
We will demonstrate this only in one dimension; however, the
generalization to 2 and 3 dimensions is obvious…
5
Method of Lines
Discretizing in Time
With finite differences, we divided time into discrete steps, as well:
u  x, t  h   u  x, t  
 2  u  x  h, t   2u  x, t   u  x  h, t  
t
h
which allowed us to find
u  x, t  h   u  x , t  
t
h
2
u  x  h, t   2u  x, t   u  x  h, t  
6
Method of Lines
Review of Finite Differences
We are given the initial state of the system by uinit(x)
Divide the space interval into n points with n – 1 intervals
7
Method of Lines
Review of Finite Differences
We evaluate the initial state at each of the n points
8
Method of Lines
Review of Finite Differences
Next, given these initial values, we take the finite-difference equation
to approximate the state at the next time t2 = t1 + t equatio
– The boundary values are defined by functions a(t) and b(t)
9
Method of Lines
Review of Finite Differences
Now, with this approximation, we approximate the values at times
t3, t4, etc.
10
Method of Lines
Method of Lines
As an alternative approach, associate with each spatial point an
unknown function uk(t)
– Two exceptions:
u1(t) = a(t)
un(t) = b(t)
This approach was popularized by the chemical engineer William E.
Schiesser in his 1991 text The Numerical Method of Lines
11
Method of Lines
Method of Lines
In order to substitute uk(t) into our mixed partial-/finite-difference
equation,


u  x, t   2  u  x  h, t   2u  x, t   u  x  h, t  
t
h
we note that the solution at location x – h is uk – 1(t) and the solution at
x + h is uk + 1(t):
d

uk  t   2  uk 1  t   2uk  t   uk 1  t  
dt
h
We also have the initial condition:
uk(tinitial) = uinit(xk)
12
Method of Lines
Systems of IVPs
Therefore, we have a system of n – 2 initial-value problems but with
n unknown functions:
d

u2  t   2  u1  t   2u2  t   u3  t  
u2  tinitial   uinit  x2 
dt
h
d

u3  t   2  u2  t   2u3  t   u4  t  
u3  tinitial   uinit  x3 
dt
h
d

u7  t   2  u6  t   2u7  t   u8  t  
dt
h
d

u8  t   2  u7  t   2u8  t   u9  t  
dt
h
u7  tinitial   uinit  x7 
u8  tinitial   uinit  x8 
13
Method of Lines
Systems of IVPs
There are two unknown functions, however, these are given by our
boundary conditions:
d

u2  t   2  u1  t   2u2  t   u3  t  
u2  tinitial   uinit  x2 
dt
h
d

u3  t   2  u2  t   2u3  t   u4  t  
u3  tinitial   uinit  x3 
dt
h
d

u7  t   2  u6  t   2u7  t   u8  t  
dt
h
d

u8  t   2  u7  t   2u8  t   u9  t  
dt
h
u7  tinitial   uinit  x7 
u8  tinitial   uinit  x8 
14
Method of Lines
Systems of IVPs
We therefore have a system of n – 2 initial-value problems
– With n = 9:
d

u2  t   2  a  t   2u2  t   u3  t  
dt
h
d

u3  t   2  u2  t   2u3  t   u4  t  
dt
h
d

u7  t   2  u6  t   2u7  t   u8  t  
dt
h
d

u8  t   2  u7  t   2u8  t   b  t  
dt
h
u2  tinitial   uinit  x2 
u3  tinitial   uinit  x3 
u7  tinitial   uinit  x7 
u8  tinitial   uinit  x8 
15
Method of Lines
Systems of IVPs
Note that we can rewrite the differential equations:
 u2  t  


 u3  t  

u t   


 u7  t  
 u t  
 8 
Thus,
 u21  t  
 1 
 u3  t  
d


u  t   u1  t   

dt
 u 1  t  
 7

1


 u t  
 8

16
Method of Lines
Systems of IVPs
We can therefore write this as:
  2

 1

1
u   t   f  t , u  t    2  
h 



where
  2
1

2
 1
def

f t, u   2  
1
h 



1
2
1
1
.2
.
1
.
1
.2
.
1
.
 a t   




0




 u t   



1 
 0 
 b t   
2 


 a t   




0



u   0 



1 
 0 
 b t   
2 


17
Method of Lines
IVP Solvers from Previous Courses
From MATH 211, you implemented the function
function [ ts, ys ] = dp45n( f, t_int, y1, h, eps_step)
where tint = [tinitial, tfinal] that uses the adaptive Dormand-Prince
method to approximate a solution to a system of initial-value
problems
y 1  t   f  t , y  t  
y  t0   y1
where we allow a maximum error of estep where we start with an
initial step size in time of h starting at tinitial and going to tfinal
18
Method of Lines
IVP Solvers from Previous Courses
We can write the function
function du = f( t, u )
M = diag( -2*ones( n - 2, 1 ) ) + ...
diag(
ones( n - 3, 1 ), 1 ) + ...
diag(
ones( n - 3, 1 ), -1 );
du = (kappa/h^2)*(M*u + [a(t); zeros( n – 4, 1 ); b(t)]);
end
  2

 1

f t, u   2  
h 



1
2
1
1
.2
.
1
.
 a t   




0



u   0 



1 
 0 
 b t   
2 

  19
Method of Lines
Example: Dirichlet Conditions
To give a specific example, consider the following heatconduction/diffusion problem:
1 u
  2u
0.2 t
that is,  = 0.2, on [0, 1] with
def
uinit  x   1
0
a t   
2
def 1

b t   
0
def
t 1
t 1
t2
t2
with n = 9 points
20
Method of Lines
Example: Dirichlet Conditions
To give this a description:
– A bar is initial uniform at 1
– For the first second, the left-hand side is placed
in contact with a heat sink at 0, after which it
is switched with a heat source of 2
– For the first two seconds, the right-hand side
is in contact with a body also at 1, after which
it is switched with a heat sink at 0
def
uinit  x   1
0
a t   
2
def 1

b t   
0
def
t 1
t 1
t2
t2
21
Method of Lines
Example: Dirichlet Conditions
Thus, we have:
function du = f111( t, u )
n = 9;
M = diag( -2*ones( n - 2, 1 ) ) + ...
diag(
ones( n - 3, 1 ), 1 ) + ...
diag(
ones( n - 3, 1 ), -1 );
kappa = 0.2;
h = 1/8;
a = @(t)(2*(t > 1));
b = @(t)(t < 2);
du = (kappa/h^2)*(M*u + [a(t); zeros( n - 4, 1 ); b(t)]);
end
22
Method of Lines
Example: Dirichlet Conditions
This can be driven by:
u_init = @(x)(x*0 + 1);
a = @(t)(2*(t > 1));
b = @(t)(t < 2);
n = 9;
xs = linspace( 0, 1, n )';
[ts, us] = dp45( @f111, [0, 3], u_init( xs(2:end - 1) ), 1, 0.0001);
us = [a(ts); us; b(ts)];
mesh( ts, xs, us );
23
Method of Lines
Example: Dirichlet Conditions
With two orientations, the output is the plot:
24
Method of Lines
Example: Dirichlet Conditions
Note that t is larger while the system converges but becomes
smaller at the discontinuities
25
Method of Lines
Example: Insulated Conditions
To give a specific example, consider the following heatconduction/diffusion problem:
1 u
  2u
0.2 t
that is,  = 0.2, on [0, 1] with
def
uinit  x   1
0 t  1
a t   
1 t  1
and an insulated right-hand boundarywith n = 9 points
def
26
Method of Lines
Example: Insulated Conditions
Using the 2nd-order divided difference approximation of the
derivative:
3u  x, t   4u  x  h, t   u  x  2h, t 

u t, x  
x
2h
and substituting our approximations for u(b, t), u(b – h, t), and
u(b – 2h, t) equate this to zero, we have:
3un  t   4un 1  t   un 2  t 
0
2h
or
4un1  t   un2  t 
un  t  
3
27
Method of Lines
Example: Insulated Conditions
We could then substitute
un  t  
4un1  t   un2  t 
3
into the last equation
d

un 1  t   2  un 2  t   2un 1  t   un  t  
dt
h
to get
4un 1  t   un  2  t  
d
 
un 1  t   2  un  2  t   2un 1  t  

dt
h 
3

 2
2

 2  un  2  t   un 1  t  
h 3
3

28
Method of Lines
Example: Insulated Conditions
Thus, we have:
function du = f112( t, u )
n = 9;
M = diag( -2*ones( n - 2, 1 ) ) + ...
diag(
ones( n - 3, 1 ), 1 ) + ...
diag(
ones( n - 3, 1 ), -1 );
M(n - 2, n - 2) = -2/3;
M(n - 2, n - 3) = 2/3;
0 t  1
Recall: a  t   
1 t  1
def
kappa = 0.2;
h = 1/8;
a = @(t)(t > 1);
du = (kappa/h^2)*(M*u + [a(t); zeros( n - 4, 1 ); 0]);
end
29
Method of Lines
Example: Insulated Conditions
This can be driven by:
u_init = @(x)(x*0 + 1);
a = @(t)(t > 1);
n = 9;
xs = linspace( 0, 1, n )';
[ts, us] = dp45( @f111, [0, 3], u_init( xs(2:end - 1) ), 1, 0.0001);
us = [a(ts); us; (4*us(end,:) - us(end - 1,:))/3];
mesh( ts, xs, us );
un  t  
4un1  t   un2  t 
3
30
Method of Lines
Example: Insulated Conditions
With two orientations, the output is the plot:
31
Method of Lines
Example: Neumann Conditions
Suppose that the right-hand side is a more general Neumann
condition:

u  t , b   ub  t 
x
Again, use the approximation
3u  x, t   4u  x  h, t   u  x  2h, t 

u t, x  
x
2h
and substituting our approximations for u(b, t), u(b – h, t), and
u(b – 2h, t) equate this to zero, we have:
3un  t   4un 1  t   un 2  t 
 ub  t 
2h
or
4u  t   un 2  t   2hub  t 
un  t   n 1
3
32
Method of Lines
The Wave Equation
We can also have an initial-value problem for the wave equation:
d2
c2
u t  2  uk 1  t   4uk  t   uk 1  t  
2 k  
dt
h
uk  tinitial   uinit  xk 
uk 1  tinitial   uinit  xk 
This can, of course, be generalized to two and three dimensions:
– There is one unknown function for each point
d

u j ,k  t   2  4u j ,k  t   u j 1,k  t   u j 1,k  t   u j ,k 1  t   u j ,k 1  t  
dt
h
 xj 
u j ,k  tinitial   uinit    
  yk  
33
Method of Lines
The Wave Equation
Recall that our function dp45.m only works with a 1st-order
differential equation
– We must therefore convert the 2nd-order ODE into a 1st-order ODE
w2 k 1  t   uk  t 
w2 k  t   uk 1  t 
34
Method of Lines
The Wave Equation
We therefore have a system of 2(n – 2) initial-value problems
– With n = 9:
d
w2  t   wn  2  t 
dt
d
w3  t   wn 3  t 
dt
d
wn  2  t   w2 n 2  t 
dt
d
wn 1  t   w2 n 1  t 
dt
d
c2
wn  2  t   2  a  t   2w2  t   w3  t  
dt
h
d
c2
wn 3  t   2  w2  t   2w3  t   w4  t  
dt
h
d
c2
w2 n  2  t   2  wn 3  t   2wn  2  t   wn 1  t  
dt
h
d
c2
w2 n 12  t   2  wn  2  t   2wn 1  t   b  t  
dt
h
w2  tinitial   uinit  x2 
w3  tinitial   uinit  x3 
wn 2  tinitial   uinit  xn 2 
wn 1  tinitial   uinit  xn 1 
wn  2  tinitial   uinit  x2 
wn 3  tinitial   uinit  x3 
w2 n 1  tinitial   uinit  xn 2 
w2 n 1  tinitial   uinit  xn 1 
35
Method of Lines
The Wave Equation
Note that we can rewrite the differential equations:
 w2  t    u2  t  

 u3  t 

w
t


 3
 


 

  u t  

 wn  2  t    n  2
 w  t    un 1  t  

   1
w  t    n 1
 wn  2  t    u2  t  
 w  t    u 1  t  

 n 3
  3


 


 
1
 w2 n  2  t    un  2  t  

 w  t   
1

 2 n 1   un 1  t  
36
Method of Lines
The Wave Equation
c2
We can therefore write this with r  2 :
h
 0








w 1  t   f  t , w  t    
2r

 r







1
0
1
0
1
0
0
r
2r
r
0
r
2r
r
0
r
0
2r
r
r
2r
 0 




0



 0 







 0 

1



0
1


w
t



 ra  t  




 0 

 0 










0
0






0
 ra  b  
37
Method of Lines
The Wave Equation
Where the function f is defined as:
 0







def

f t, w   
2r

 r







1
0
1
0
1
0
0
r
2r
r
0
r
2r
r
0
r
0
2r
r
r
2r
 ra  t  




0



 0 







 0 

1



1
 ra  b  
w

 0 




 0 




0









0
0






0
 0 
38
Method of Lines
The Wave Equation
Thus, we have:
function [dw] = f113( t, w )
n = 17;
c = 1;
h = 1/8;
r = (c/h)^2;
Id = ones( n - 2, 1 );
Z = zeros( n - 2, n - 2 );
Id_off = ones( n - 3, 1 );
M = diag( -2*r*Id ) + diag( r*Id_off,
M = [Z eye( n - 2, n - 2 ); M Z];
1 ) + diag( r*Id_off, -1 );
v = zeros( 2*n - 4, 1 );
a = @(t)(sin(3*t));
b = @(t)(0*t);
v(n - 1) = r*a(t);
v(end) = r*b(t);
dw = M*w + v;
end
39
Method of Lines
The Wave Equation
This can be driven by:
u_init = @(x)(x*0);
du_init = @(x)(x*0);
a = @(t)(sin(3*t));
b = @(t)(0*t);
n = 17;
xs = linspace( 0, 1, n )';
U_init = [u_init( xs(2:end - 1) ); du_init( xs(2:end - 1) )];
[ts, us] = dp45( @f113, 0, 10, U_init, 1, 0.0001);
size( us )
u = [a(ts); us(1:(n - 2), :); b(ts)];
mesh( ts, xs, u );
du = [cos(3*ts); us((n - 1):end, :); 0*ts];
% Plot the derivative...
mesh( ts, xs, du );
40
Method of Lines
The Wave Equation
The plot of the solution and the derivative are here:
41
Method of Lines
Finite Element Methods
It will also work equally effectively for finite element methods:
– The elements define the relationship between a function defined at a
point and its neighbouring functions
42
Method of Lines
Summary
In this topic, we have covered the method of lines:
– We discretize only the space component
– This creates a system of initial-value problems
– This may be solved using a high-accuracy adaptive ODE solver
43
Method of Lines
References
[1] Schiesser, W. E., The Numerical Method of Lines, Academic
Press, 1991.
44
Method of Lines
Usage Notes
•
•
These slides are made publicly available on the web for anyone to use
If you choose to use them, or a part thereof, for a course at another
institution, I ask only three things:
– that you inform me that you are using the slides,
– that you acknowledge my work, and
– that you alert me of any mistakes which I made or changes which you make, and
allow me the option of incorporating such changes (with an acknowledgment) in
my set of slides
Sincerely,
Douglas Wilhelm Harder, MMath
[email protected]
45