ENGR90024 COMPUTATIONAL FLUID DYNAMICS
Lecture O02 Euler’ss method, programming and the Unix environment Euler’
Lecture O01
φ
= 1
− e−t
1
0.9
0.8
dφ = f (t, φ) dt tmin
s h o d t h e M ic a l t i c a h e m t h M a
Initiall Value Problem Initia Pr oblem
0.6
) t (
!
0.5
0.4
0.3
0.2
0.1
≤ t ≤ tmax
φ(tmin ) = φmin
0.7
Exact solution
Compare N u m e r i ic c a l / /C o m p u t e er M e t h ho d s
Started in Lecture O01 will finish it in Lecture O02
Approximate solution
0 0
1
2
3
4 t
5
6
7
8
dφ = f (t, φ) dt tmin
≤ t ≤ tmax
φ(tmin ) = φmin
Initial Value Problem
N u m e r i ic c a l / /C o m p u t e er M e t h ho d s
Started in Lecture O01 will finish it in Lecture O02
Approximate solution
dφ = f (t, φ) dt min
≤
≤
max
φ(tmin) = φ min
N u m e r i ic c a l / / C o m p u t e e r M e t h h o d s
Initial Value Problem
Approximate solution
Explicit Euler Method (section 1.2.1 of the printed lecture notes) To derive Euler’s method, start with Eq. (O01.7) l+1
φ(t
l
) = φ(t ) +
∆t dφ
1!
dt
l
(t )
If we truncate tr uncate Taylor aylorss series, ser ies, Eq. (O01.7) (O01.7),, and ignore term te rmss of the order of ! t2, we will get the following approximation for the value of "(tl+1) (remember that tl+1=tl+!), l+1
φ(t
l
l
) = φ(t ) + ∆t l
dφ dt
l
(t )
φ(t + ∆t) = φ(t ) + ∆t
dφ dt
l
(t )
(O02.1)
Note that the above expression has error O(! t2). Recall that ODEs are typically expressed as
dφ (t) = f (φ, t) dt
(O02.2)
Substituting Eq. (O02.2) into Eq. (O02.1) l
l
φ(t + ∆t) = φ(t ) + ∆t f (φ, t) to obtain
φ
l+1
l
l
l
= φ + ∆tf (φ , t )
(O02.3)
where "l+1 is the approximate value of "(tl+1) and "l is the approximate value of "(tl).The formula above is called the explicit Euler’ss method and it can be used to approximate the solution at time Euler’ level tl+1 when you know the (approximate) solution at time level tl.
φ (t )
Smaller ! t will lead to smaller error Predicted value of "l+1 Error
dφ = f (φl , tl ) dt
φ(tl+1 )
φ(t)
True value of "l+1 φl+1 = φl + ∆tf (φl , tl )
φ(tl )
Euler’s formula
∆t l
t
l+1
t
φ (t )
Smaller ! t will lead to smaller error
Predicted value of "l+1
dφ = f (φl , tl ) dt
Error φ(t)
φ(tl+1 )
True value of "l+1 φl+1 = φl + ∆tf (φl , tl )
φ(tl )
Euler’s formula
∆t l
t
l+1
t
φ (t )
Smaller ! t will lead to smaller error
dφ = f (φl , tl ) dt
Predicted value of "l+1 φ(t) Error True value of "l+1
φ(tl+1 )
φl+1 = φl + ∆tf (φl , tl )
φl
Euler’s formula
∆t l
t
t
l+1
Predicted "(t)
φ(t)
∆t
∆t
∆t
∆t
(t3, "3) dφ(t) dt
dφ(t
∆t
dt
∆t
(t2, "2) dφ(t dt
∆t
(t4, "4) dφ(t)
φmin
dt
(t1,"1)
(t5, "5)
2
5
∆t
Example O02.1 (similar to Example 1.1 of the printed lecture notes): Using Euler’s method, solve
dφ dt
= 1
For 0 < t < 8 with "(t=0)=0 and a) ! t=2 b) ! t=1 c) ! t=0.5 d) ! t=0.1
−φ
φ(t)
∆t
∆t
∆t
t3, dφ(t) dt
dt
function MPO02p1() clear all; close all; Delta_t=2.0; t=0:Delta_t:8.0 phi(1)=0.0; for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1phi(n)) end plot(t,phi,'ko-') hold on ezplot(@(t)1-exp(-t),[0,8,0,2]) xlabel('t'); ylabel('\phi'); legend('Euler','True');
3
dφ t
∆t
dt
t2 dφ(t)
∆t
∆t
2
∆t
t4
4
dφ(t)
φmin
t1
dt
1
t5
t1=tmin
t
t
t
∆t
5
t5=tmax
function MPO02p1() clear all; close all; Delta_t=2.0; t=0:Delta_t:8.0 phi(1)=0.0;
Output
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1-phi(n)) end plot(t,phi,'ko-') hold on ezplot(@(t)1-exp(-t),[0,8,0,2]) xlabel('t'); ylabel('\phi'); legend('Euler','True');
1!exp(!t)
2
Euler True
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0
1
2
3
4 t
5
6
7
8
!
2
2
! Euler True
Euler True
1.8
1.8
!t=2.0
1.6
1.4
1.4
1.2
1.2
1
!
!t=1.0
1.6
!
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
1
2
3
4
5
6
7
0
8
t
2
0
1
2
3
2
Euler
4 t
5
6
7
Euler
True
1.8
True
1.8
!t=0.5
1.6
1.4
8
!t=0.1
1.6
1.4 1.2
1.2 !
1 !
1
0.8
0.8 0.6
0.6 0.4
0.4 0.2
0.2 0 0
1
2
3
4 t
5
6
7
8
0 0
1
2
3
4 t
5
6
7
8
End of Example O02.1
φ
= 1
− e−t
1
0.9
0.8
dφ = f (t, φ) dt tmin
d s t h o e l M t ic a a m t h e M a
Initial Value Problem
0.6
) t (
!
0.5
0.4
0.3
0.2
0.1
≤ t ≤ tmax
φ(tmin ) = φmin
0.7
Exact solution
0
Compare N u m e r ic a l /C o m p u t er M e t ho d s
0
1
2
3
4
5
6
7
8
t
2 Euler True
1.8
Approximate solution
1.6
1.4
1.2
!
1
0.8
0.6
0.4
0.2
0 0
1
2
3
4 t
5
6
7
8
Local & Global Truncation Error (page 6 printed notes) φ(t
l+1
l
) = φ( t ) +
∆t dφ
1!
dt
(tl )
Euler method, Eq. (O02.3) has local truncation error O(! t2). This means that if you halve ! t, you can expect the error to reduce by a quarter. However, we have to take O(! t-1) steps to find the solution at time tNt. A conservative position would to assume that the error accumulate at every time step E global t
N N = φ(t ) − φ 0 t
=
t
0
0
1
1
N t −1
φ(t ) − φ0 + ∆tf (φ , t ) + ∆tf (φ , t ) + . . . + + ∆tf (φ
N t −1
,t
)
Hence, the global truncation error is O(! t). In general, if a method has local truncation error of O(! tN), then the global truncation error is O(! tN-1).
This is a very inefficient MATLAB program!! function MPO02p1() clear all; close all; Delta_t=2.0; t=0:Delta_t:8.0 phi(1)=0.0; for n=1:length(t)-1 n=1:length(t)-1 for phi(n+1)=phi(n)+Delta_t*(1-phi(n)) phi(n+1)=phi(n)+Delta_t*(1-phi(n)) end end plot(t,phi,'ko-') hold on ezplot(@(t)1-exp(-t),[0,8,0,2]) xlabel('t'); ylabel('\phi'); legend('Euler','True');
How to make the our Matlab program run faster: Memory Allocation (see Example 1.1 page 7 of printed lecture notes) In Matlab program MPO02p1.m, a new block of memory is allocated to the variable phi() at every single iteration of the “for” loop for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1-phi(n)) end
This is very inefficient because Matlab has to find a new block of memory for phi() every time it goes into the for loop. It then has to copy the data to this new block. This is very inefficient. It is much more efficient to preallocate memory to phi() before you use it in the “for” loop. You can preallocate memory and zero every element in the vector by issuing the command phi=zeros(1,100);
This command allocates memory to 100 elements in the phi() vector (1 x 100) and all elements have a zero value.
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0; for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
0
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
0 "(2)
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
0
"(2)
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
0
"(2)
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
0
"(2) "(3)
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
0
"(2) "(3)
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
0
"(2) "(3)
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0; for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
0
phi=
0 "(2) "(3) "(4)
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
"(2) "(3)
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
0
"(2) "(3) "(4)
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
0
"(2) "(3) "(4)
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0; for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
0
!(2) !(3) !(4)
phi=
0
!(2) !(3) !(4) !(5)
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
0
"(2) "(3) "(4) "(5)
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
0
"(2) "(3) "(4) "(5)
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
0
0
0
0
0
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
0
"(2) "(3) "(4) "(5)
0
"(2) 0
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
0
0
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
0
"(2) "(3) "(4) "(5)
0
"(2) "(3) 0
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
0
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
0
"(2) "(3) "(4) "(5)
0
"(2) "(3) "(4) 0
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
0
"(2) "(3) "(4) "(5)
0
"(2) "(3) "(4) "(5)
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
Without preallocation of memory t=0:Delta_t:8.0 phi(1)=0.0;
phi=
0
"(2) "(3) "(4) "(5)
0
"(2) "(3) "(4) "(5)
for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
With preallocation of memory t=0:Delta_t:8.0 phi=zeros(1,length(t)); for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end
phi=
CASE STUDY O02.2 (similar to Example 1.1 of the printed lecture notes): Execute MPO02p1.m for (I) ! t=1.0e-4 (II)! t=1.0e-5 (III)! t=1.0e-6 (IV)! t=1.0e-7 For each value of ! t, and how much time it takes to run. Preallocate memory to phi() and run the same program again. See how long it now takes to execute. Use the tic() and toc() function to record the time it takes to execute a certain command in Matlab.
function MPO02p2() clear all; close all; Delta_t=1.0e-4; %Delta_t=1.0e-5; %Delta_t=1.0e-6; %Delta_t=1.0e-7; t=0:Delta_t:10.0; % preallocating memory & set elements to zero %phi=zeros(1,numel(t)); phi(1)=0.0; tic % start timer for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1-phi(n)); end temp=toc; %end timer fprintf(1,'It takes %e seconds to execute the Matlab program\n',temp);
Time it takes to execute MPO02p2.m with and without preallocation of the array phi() for different !t
End of Case Study O02.2
Unix Programming Environment
Example O02.3: ENGR90024 Write a C++ program that uses Euler’s method to solve
COMPUTATIONAL FLUID DYNAMICS
dφ
= 1
−φ
Lecture dt
O04
For 0 < t < 8 with "(t=0)=0 and Ordinary Differential Equations (ODEs): a) !t =2 b) !t =1 c) !t =0.5 C Programming & the Unix Environment d) !t =0.1
In the Unix environment, most tasks are conducted from the command line.
• Use the mkdir command to create a directory, O02p3,
for your
files
mkdir O02p3
• Check that the directory is there by using the ls command. ls O02p3
• Use the cd command to change directory to O02p3 cd O02p3
•
Edit a file called CPO02p3.cpp with the gedit program by typing
gedit CPO02p3.cpp •
Type in the following C++ program (looks very similar to the Matlab program MPO02p1.m)
C++ Program CPO02p3.cpp #include #include #include #include #include
using namespace std; int {
main(int argc, char **argv) /* Simulation parameters */ float t_min = 0.0; float t_max = 10.0; float Delta_t = float phi_min = int N_t int l fstream OutputFile;
0.5; 0.0; = (t_max-t_min)/Delta_t; = 0;
/* Allocate arrays */ float t[N_t+1]; float phi[N_t+1]; float phi_analytical[N_t+1]; /* NOTE: Here we are using 'static allocation' to store the arrays. */ /* Set initial conditions */ t[0] = t_min; phi[0] = phi_min; phi_analytical[0] = phi_min; /* Main time marching loop for(l=0; l
*/ = t[l] + Delta_t; = phi[l] + Delta_t*(1 - phi[l]); = 1 - exp(-t[l+1]);
/* Write the simulation data to an output file */ OutputFile.open("CPO02p3.out", ios::out); for(l=0; l<=N_t; l++) { OutputFile << t[l] << "\t" << phi[l] << "\t" << phi_analytical[l] << endl; } OutputFile.close(); return 0;
C++ does not have plotting capability. You need to write the data to a file and use MATLAB to plot
MATLAB Program function MPO02p1() clear all; close all; Delta_t=0.1; t=0:Delta_t:8.0 phi(1)=0.0; for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end plot(t,phi,'ko-') hold on ezplot(@(t)1-exp(-t),[0,8,0,2]) xlabel('t'); ylabel('\phi'); legend('Euler','True');
C++ Program #include #include #include #include #include using namespace std; int main(int argc, char **argv) { /* Simulation parameters */ float t_min = 0.0; float t_max = 10.0; float Delta_t = float phi_min = int N_t int l fstream OutputFile;
0.5; 0.0; = (t_max-t_min)/Delta_t; = 0;
/* Allocate arrays */ float t[N_t+1]; float phi[N_t+1]; float phi_analytical[N_t+1]; /* NOTE: Here we are using 'static allocation' to store the arrays. */ /* Set initial conditions */ t[0] = t_min; phi[0] = phi_min; phi_analytical[0] = phi_min; /* Main time marching loop for(l=0; l
*/ = t[l] + Delta_t; = phi[l] + Delta_t*(1 - phi[l]); = 1 - exp(-t[l+1]);
/* Write the simulation data to an output file */ OutputFile.open("CPO02p3.out", ios::out); for(l=0; l<=N_t; l++) { OutputFile << t[l] << "\t" << phi[l] << "\t" << phi_analytical[l] << endl; } OutputFile.close(); return 0;
MATLAB Program function MPO02p1() clear all; close all; Delta_t=0.1; t=0:Delta_t:8.0 phi(1)=0.0; for n=1:length(t)-1 phi(n+1)=phi(n)+Delta_t*(1 phi(n)) end plot(t,phi,'ko-') hold on ezplot(@(t)1-exp(-t),[0,8,0,2]) xlabel('t'); ylabel('\phi'); legend('Euler','True');
C++ Program #include #include #include #include #include using namespace std; int main(int argc, char **argv) { /* Simulation parameters */ float t_min = 0.0; float t_max = 10.0; float Delta_t = float phi_min = int N_t int l fstream OutputFile;
0.5; 0.0; = (t_max-t_min)/Delta_t; = 0;
/* Allocate arrays */ float t[N_t+1]; float phi[N_t+1]; float phi_analytical[N_t+1]; /* NOTE: Here we are using 'static allocation' to store the arrays. */ /* Set initial conditions */ t[0] = t_min; phi[0] = phi_min; phi_analytical[0] = phi_min; /* Main time marching loop for(l=0; l
*/ = t[l] + Delta_t; = phi[l] + Delta_t*(1 - phi[l]); = 1 - exp(-t[l+1]);
/* Write the simulation data to an output file */ OutputFile.open("CPO02p3.out", ios::out); for(l=0; l<=N_t; l++) { OutputFile << t[l] << "\t" << phi[l] << "\t" << phi_analytical[l] << endl; } OutputFile.close(); return 0;
•
“compile” the C++ program by typing the command
g++ -o CPO02p3 CPO02p3.cpp
• Run/execute the program by typing ./CPO02p3
• After you execute the program, you should find that you have the file
CPO02p3.out in that directory
•
Check this by typing the ls command
ls CPO02p3 CPO02p3.cpp CPO02p3.out
• Execute the Matlab program “CPO02p3_Postprocessing.m” to plot the the results (see next slide).
%%%%%%%%%%%%%%%%%%%%%%%%%%% % CPO02p3_Postprocessing.m %%%%%%%%%%%%%%%%%%%%%%%%%%% close all, clear all, clc; load 'CPO02p3.out' [N_t N_e]
= size(CPO02p3);
% Plot the solution plot(CPO02p3(:,1), CPO02p3(:,2), '.-', CPO02p3(:,1), CPO02p3(:,3), '-', 'LineWidth', 2, 'MarkerSize', 20); axis([0 10 0 2]); grid on; xlabel('t'); ylabel('\phi');
2
1.8
1.6
1.4
1.2
!
1
0.8
0.6
0.4
0.2
0 0
1
2
3
4
5 t
6
7
8
9
10
Flow diagram
CPO02p3.cpp g++ CPO02p3
CPO02p3.out CPO02p3_Postprocessing.m 2
1.8
1.6
1.4
1.2
!
1
0.8
0.6
0.4
0.2
