Numerical Methods in Engineering (ENGR20005)
Lecture 17 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L16.3:
Solving Boundary Value Problem Using Spectral Methods
2
Lecture Recap
y0 = α
D2 D2 … D2
11 12 1,n−1
D2 D2 … D2
21 22 2,n−1
⋮⋮…⋮
D2 D2 D2 D2
n−1,1 n−1,2 43 n−1,n−1
e4x1 e4x2
⋮ e4xn−1
D2 1n
β D2 2n
⋮ D2 n−1,n
Lecture Recap
y1 y2 yn−2 yn−1
yn = β x0=−1 x1 x2 xn−2 xn−1 xn=1
So in general if you split your domain up into n intervals, you will need to solve an equation that looks like
D2 y1 10
y2 D2 ⋮= −α20−
y⋮ n−1 D2
Two Dirichlet Boundary Condition
3
n−1,0
y0 = α
Lecture Recap
y1 y2 yn−2 yn−1
x0=−1 x1 x2 xn−2 xn−1 xn=1
y′ = β n
So in general if you split your domain up into n intervals, you will need to solve an equation that looks like
D2 D2 D2 11 12 13
D2 D2 D2 21 22 23
D2 D2 D2 31 32 33
D2 D2 D2 41 42 43
Dn,1 Dn,2 Dn,3
… D2 D2 1,n y1 e4x1 10
…D2 2 2,n y2 e4x2 D20
… ⋮ ⋮ = ⋮ −α ⋮ yn−1 e4xn−1 D2
… D2 y n−1,0 n−1,n n β D
… Dn,n n0
One Dirichlet, one Neumann Boundary Condition
4
What if the interval is NOT from -1 to 1???
L17.1:
Boundary Value Problem: Spectral Method with GLo points in
domain other than −1 ≤ x ≤ 1
Example L17.1: Solve the following linear BVP using spectral differencing schemes. Use the Lagrange basis functions and the GLo points. Taken from “Applied Numerical Methods Using MATLAB” by Yang, Cao, Chung and Morris
y′′(t) + 2t y′(t) − t2 y(t) = 0 2
t ∈ [1,2] and the use boundary conditions y(1) = 5 and y(2) = 3.
y(t) = t + t4 2
Compare your solution with the exact expression for
This boundary problem is expressed in the domain 1 ≤ t ≤ 2 but the GLo points are usually only defined −1 ≤ x ≤ 1. So we need to transform the problem above from t ∈ [1,2] to x ∈ [−1,1]. In
lecture 16, it was shown that to transform from t ∈ [a, b] to x ∈ [−1,1] we need to make the change of
variables
a=1
t = b − a x + b + a b = 2 ⇒ t = 1 (x + 3)
222
dy=( 2 )dy dt b−a dx
d2y=( 2 )2d2y dt2 b−a dx2
6
⇒ dy = 2 dy dt dx
⇒ d2y = 4 d2y dt2 dx2
y′′(t) + 2t y′(t) − t2 y(t) = 0 t ∈ [1,2] ⇒ t = 12 (x + 3) 2⇒
4y′′(x)+ 8 y′(x)− 8 y(x)=0 x+3 (x+3)2
⇒ y′′(x)+ 2 y′(x)− 2 y(x)=0 x∈[−1,1]
x+3 (x+3)2 The boundary conditions
y(t=1)=5⇒ y(x=−1)=5 y(t=2)=3⇒ y(x=1)=3
To find the structure of the matrix system we need to solve, we will look at the domain −1 ≤ x ≤ 1 discretise into 5 intervals.
7
dy =2dy dt dx
d2y =4d2y dt2 dx2
y(x = − 1) = 5
y0 = α y1
x0=−1 x1
y′′
y2
x2 2
y3 y4 2×3 x4
y(x = 1) = 3 y5 = β
x5=1
y′′(x) + x + 3 y′(x) − (x + 3)2 y(x) = 0 (2/(x +3))y′ (2/(x +3)2)y
000000
y′′ (2/(x +3))y′ (2/(x +3)2)y 111110
y′′ (2/(x +3))y′ (2/(x +3)2)y
2+ 2 2− 2 2=0
y′′ (2/(x +3))y′ (2/(x +3)2)y 333330
y′′ (2/(x +3))y′ (2/(x +3)2)y 444440
y′′ (2/(x +3))y′ (2/(x +3)2)y 55555
8
y′′ 0
y′′ 1
y′′ 2
y′′ 3
y′′ 4
y′′ 5
( 2 / ( x 0 + 3 ) ) y 0′ ( 2 / ( x 1 + 3 ) ) y 1′ ( 2 / ( x 2 + 3 ) ) y 2′ ( 2 / ( x 3 + 3 ) ) y 3′ ( 2 / ( x 4 + 3 ) ) y 4′ ( 2 / ( x 5 + 3 ) ) y 5′
(2/(x0 + 3)2)y0 (2/(x1 + 3)2)y1 (2/(x2 + 3)2)y2 (2/(x3 + 3)2)y3 (2/(x4 + 3)2)y4 (2/(x5 + 3)2)y5
{y′′} = [D2]{y}
y′′ 0
y′′ 1
y′′ 2
y′′ 3
y′′ 4
y′′ 5
D2D2D2D2D2D2 000102030405 y0
D2D2D2D2D2D2 y 101112131415 1
+−=
0
{(2/(x + 3))y′} = diag((2/(x + 3)))[D]{y}
{(2/(x + 3)2)y} = diag((2/(x + 3)2)){y}
[[D2] + diag {2/(x + 3)}[D] − diag {2/(x + 3)2}] {y} = {0}
D2 D2 D2 D2 D2 D2 y =202122232425 2
D2 D2 D2 D2 D2 D2 y3 303132333435 y4
D2D2D2D2D2D2 y 404142434445 5
D2D2D2D2D2D2 505152535455
(2/(x0 + 3))y0′ 2/(x0 + 3)
(2/(x1+3))y1′ 0
(2/(x + 3))y′ 0 22=
(2/(x3+3))y3′ 0 (2/(x4+3))y4′ 0 (2/(x5+3))y5′ 0
0 2/(x1 + 3) 0
0
0
2/(x +3) 2
0 0 0
0 0 0
0
0
0
0
0 2/(x5 + 3)
D00 D01 D D
D02 D03 D D2
D04 D05 y0 D D y
where
[A]{y} = {0}
[A] = [D2] + diag {2/(x + 3)}[D] − diag {2/(x + 3)2}
0 0 2/(x3+3) 0 000 2/(x4 + 3) 000 0
D10 D11 20 21 D30 D31 D40 D41 D50 D51
D12 D13 22 23 D32 D33 D42 D43 D52 D53
D14 D15 y1 24 25 y2
D34 D35 y3 D44 D45 y4 D54 D55 5
(2/(x0 + 3)2)y0 (2/(x1 + 3)2)y1 (2/(x2 + 3)2)y2 (2/(x3 + 3)2)y3 (2/(x4 + 3)2)y4 (2/(x5 + 3)2)y5
=
2/(x0+3)2 0 0 0 0 0
0 2/(x1+3)2 0 0 0 0
0 0 2/(x2+3)2 0 0 0
0 0 0 2/(x3+3)2 0 0
0 0 0 0 2/(x4+3)2 0
0 0 0 0 0 2/(x5+3)2
y0 y1 y2 y3 y4 y5
9
y(x = − 1) = 5
y0 = α y1
y2 x2
y3 y4
y(x = 1) = 3 y5 = β
x0=−1 x1
y′′(x)+ 2 y′(x)− 2 y(x)=0
x3
x4 x5=1 [A]{y} = {0}
x+3 (x+3)2
y′′ (2/(x +3))y′ (2/(x +3)2)y
000000
y′′ (2/(x +3))y′ (2/(x +3)2)y 111110
y′′ (2/(x +3))y′ (2/(x +3)2)y
2+ 2 2− 2 2=0
y′′ (2/(x +3))y′ (2/(x +3)2)y 333330
y′′ (2/(x +3))y′ (2/(x +3)2)y 444440
y′′ (2/(x +3))y′ (2/(x +3)2)y 55555
y0 y1
y2 = y3
y4 y5
0 0 0 0 0 0 10
[A] = [D2] + diag {2/(x + 3)}[D] − diag {2/(x + 3)2}
A00 A01 A02 A03 A04 A05 A10 A11 A12 A13 A14 A15 A20 A21 A22 A23 A24 A25 A30 A31 A32 A33 A34 A35 A40 A41 A42 A43 A44 A45 A50 A51 A52 A53 A54 A55
y(x = − 1) = 5 y0 = 5
y1y2 y3y4
y(x =1)=3 y5 = 3
0 0
100000
A10 A11 A12 A13 A14 A15 A20 A21 A22 A23 A24 A25 A30 A31 A32 A33 A34 A35 A40 A41 A42 A43 A44 A45
000001
x0=−1 x1
x1:
x2: x3: x4: x5:
A10 A11 A12 A20 A21 A22 A30 A31 A32 A40 A41 A42 A50 A51 A52
A03 A04 A05
A13 A14 A15 y1
y′′(x)+ 2 y′(x)− 2 y(x)=0
x2
x0: A00 A01 A02
x3 x4 x5=1 y0 0
A23 A24 A25
A33 A34 A35 y3
x+3 (x+3)2
y2 = A43 A44 A45 y4
A53 A54 A55
y5 0 y0 5
11
y1
y2 =
y3
y4 53
y0 = 5
y1y2 y3y4
y5 = 3
0 0 0 0
x0=−1 x1 x2
x3 x4 x5=1 0
A00 A
A01 A02 A A
A03 A04 A05
y 0
10 A
11 12 A A
A A A
13 14 15 y
20 A
21 22 A A
1 A A A y
30 A
31 32 A A
33 34 35 A A A
4
40 A50
41 42 A51 A52
43 44 45 A53 A54 A55
y 5
0
23 24 25 2
y A A A 3
= y
12
y0 = 5
y1y2 y3y4
y5 = 3
0
x0=−1 x1
x2 x3 x4 x5=1
A00 A01 A02 A03 A04 A05 y0 0 A10 A11 A12 A13 A14 A15 y1
A20 A21 A22 A23 A24 A25 y2 =
A30 A31 A32 A33 A34 A35 y3
A40 A41 A42 A43 A44 A45 y4
A50 A51 A52 A53 A54 A55 y5 0
AAAAAA y0
10 11 12 13 14 15 y1 0
A20 A21 A22 A23 A24 A25 y2 =0 A30 A31 A32 A33 A34 A35 y3 0
A40 A41 A42 A43 A44 A45 y4 0 y5
13
y0 = 5
y1y2 y3y4
y5 = 3
A10 A20 A30 A40
A15 A25 A35 A45
A10 A20 A30 A40
A15 A25 A35
A45
x0=−1 x1
x2
A11 A12 A13 A14
x3 =0
x4 x5=1
y y0
0
AAAA
y1
21 22 23 24 A31 A32 A33 A34 A41 A42 A43 A44
y2
0 0
A11 A12 A13 A14
y1
y3 4
y5
0 y +A21A22A23A24y2+y =0 0 A31A32A33A34 y3 5 0 y4 0
A41 A42 A43 A44
14
y0 = 5
y1y2 y3y4
y5 = 3
A10 A20 A30 A40
A15 A25 A35
A45
A10 A20 A30 A40
A15 A25 A35
A45
x0=−1 x1 x2 x3 x4 x5=1 A11 A12 A13 A14 y1 0
y +A21A22A23A24y2+y =0
0 A31A32A33A34 y3 5 A41 A42 A43 A44 y4
0 0
Rearranging gives
A11 A12 A13 A14 y1 A21 A22 A23 A24 y2 A31 A32 A33 A34 y3 A41 A42 A43 A44 y4
0
= 0 0
0
−5
−3
Can solve to find y1, y2, y3 and y4. 15
y0 = α
y1 y2 yn−2 yn−1
yn = β
0 ⋮0
A10 A20
A11 A12 … A1,n−1 A21 A22 … A2,n−1
⋮⋮…⋮ An−1,1 An−1,2 A43 An−1,n−1
α
n−1,0
A⋮
β
⋮ An−1,n
A1n
A2n
x0=−1 x1 x2 xn−2 xn−1 xn=1
So in general if you split your domain up into n intervals, you will need to solve an equation that looks like
y1 y2
=−− n−1
y⋮
16
Lecture17M01.m
clear all y0′
Npoints=50; y1′
y0 y1
alpha = 5; %BC1 y(x=-1)=5 y2′ = y2
beta = 3; %BC2 y(x=1)=3 ⋮ y⋮
y n′ − 1 %xpoints=linspace(-1,1,Npoints)’; yn′
n − 1 yn
[xpoints,w,P]=lglnodes(Npoints); xpoints=-xpoints;
1 y′′
tpoints=(xpoints+3)/2;
D=DerivMatrix(xpoints,n);
D2 = D*D;
A=D2+diag(2./(xpoints+3))*D-diag(2./(xpoints+3).^2);
Amatrix=A(2:end-1,2:end-1);
Cvector=-alpha*A(2:end-1,1)-beta*A(2:end-1,end);
t = 2 (x + 3) 0 y0 y′′
y = @(t) t+4./t.^2; yexact = y(tpoints);
n=length(xpoints)-1; %order of polynomial
1 y1 y′′ y2
⋮2= ⋮
y′′ n−1
y′′ n
yn−1 yn
A11 A12 … A1,n−1 A21 A22 … A2,n−1
⋮⋮…⋮ An−1,1 An−1,2 A43 An−1,n−1
y0 yα y1 y1
y2 = ⋮2 y⋮ yn−1
n−1
yn β
yapprox = [alpha;sol;beta];
D00 D01 D02 D03 … D0,n
D10 D11 D12 D13 … D1,n
D20 D21 D22 D23 … D2,n ⋮⋮⋮⋮…⋮
Dn−1,0 Dn−1,1 Dn−1,2 Dn−1,3 . . . Dn−1,n Dn,0 Dn,1 Dn,2 Dn,3 . . . Dn,n
D2 D2 D2 D2 … D2 00 01 02 03 0,n
D2 D2 D2 D2 … D2 10 11 12 13 1,n
D2 D2 D2 D2 … D2 20 21 22 23 2,n
⋮⋮⋮⋮…⋮
D2 n−1,0
D2 n,0
D2 D2 D2 n−1,1 n−1,2 n−1,3
D2 D2 D2 n,1 n,2 n,3
… D2 n−1,n
… D2 n,n
[A] = [D2] + diag {2/(x + 3)}[D] − diag {2/(x + 3)2}
A10
α A20 −β
A1n
y1
n−1 0
sol = Amatrix\Cvector;
hold off plot(tpoints,yapprox,’bo’) hold on plot(tpoints,yexact,’k-‘) xlabel(‘t’);ylabel(‘y’);
function D=DerivMatrix(x,n) D=zeros(n+1,n+1);
for i=1:n+1
0
y⋮2 = ⋮0 −
A⋮ n−1,0
A2n
⋮ An−1,n
17
num(i)=1.0;
A11 A12 … A1,n−1 A21 A22 … A2,n−1
⋮⋮…⋮ An−1,1 An−1,2 A43 An−1,n−1
A10
α A20 −β
A1n
y1 0
yapprox = [alpha;sol;beta];
y2 = 0 − y⋮ ⋮
A⋮ n−1,0
A2n
⋮ An−1,n
n−1 0
hold off plot(tpoints,yapprox,’bo’) hold on plot(tpoints,yexact,’k-‘) xlabel(‘t’);ylabel(‘y’);
y0 yα y1 y1
y2 = ⋮2 y⋮ yn−1
n−1
yn β
Approximated solution
Exact solution y(t) = t + t4 2
sol = Amatrix\Cvector;
function D=DerivMatrix(x,n) D=zeros(n+1,n+1);
for i=1:n+1
num(i)=1.0; for k=1:n+1
if k ~= i num(i)=num(i)*(x(i)-x(k));
end end
end
for j=1:n+1 den(j)=1.0; for k=1:n+1
if k ~= j den(j)=den(j)*(x(j)-x(k));
end end
end
for i=1:n+1 for j=1:n+1
if i ~= j D(i,j)=num(i)/(den(j)*(x(i)-x(j)));
end end
end
for i=1:n+1 for k=1:n+1
End of Example L17.1
L17.2:
Ordinary Differential Equations: Initial Value Problem
20
Some examples of initial value problems.
l
θ
The motion of the ball is governed by
d2θ + g sin θ = 0 dt2 l
21
Ordinary Differential Equations
In many engineering problems, you will need to solve differential equations that look something like
(17.1)
dx dt
in the domain
with the initial condition
= f(t, x) a≤t≤b
x(t = a) = α
This is called the initial value problem where all the information you are
given about x is at one end (usually at t = a) of the domain. 22
x
Solve
dx dt
= f(t, x)
to find value of x(t) for
dx dt
= f(t, x) a≤t≤b
x(t = a) = α
a≤t≤b
Initial Value Problem
α
a
t
b
Exact (true) solution
dx dt
= f(t, x) a≤t≤b
x(t = a) = α
Compare
Initial Value Problem
Approximate solution
Mathematical Methods
Numerical/Computer Methods
Exact (true) solution
dx dt
= f(t, x) a≤t≤b
x(t = a) = α
Compare
Approximate solution
Initial Value Problem
Example L17.2:
Solve the following ordinary differential equation (ODE)
dx dt
=1−x
(17.2)
with x(t = 0) = 0 Use MATLAB to plot the solution for 0 ≤ t ≤ 8.
Mathematical Methods
Numerical/Computer Methods
=1−x Rearranging gives
dx dt
dx dt
+x=1
Use integrating factor, multiply both sides by et
tdx t t e dt + e x = e 1
d
(etx) = et Integrating both sides with respect to t gives
dt
d dt
(etx) = et Integrating both sides with respect to t gives
etx = et + K
where K is a constant. We are given that x = 0 when t = 0. Substituting
into the above equation gives
e0 × 0 = e0 + K 1×0=1+K
K=−1
Thus the solution to Eq. (17.2) is
x(t) = 1 − e−t
t=[0.0 0.1 0.2 0.3 0.4……..8.0]
Lecture17M02.m
t=0:0.1:8
x=1-exp(-t)
plot(t,x,’k-‘,’linewidth’,3) xlabel(‘t’)
ylabel(‘x(t)’)
1-e-0.0
x=[0.0
1-e-0.2
0.0952 0.1813
1-e-0.1
1-e-8.0
0.2592 …….. 1.0]
1-e-0.3
End of Example L17.2
x = 1 − e−t
Exact solution
dx dt
= f(t, x) a≤t≤b
x(t = a) = α
Compare
Initial Value Problem
Approximate solution
Mathematical Methods
Numerical/Computer Methods
Taylor series
Recall that the Taylor series expansion of a function x(t) about the point tn can be written as
∞ (tn+1−tn)kdkx x(tn+1) = ∑k=0 k! dtk (tn)
tn+1 = tn + Δt.
Equation (17.3) can be expanded to look like
x(tn+1)=x(tn)+(tn+1−tn)dx(tn)+(tn+1−tn)2 d2x(tn)+(tn+1−tn)3 d3x(tn)+….. 1!dt 2!dt2 3!dt3
x(tn+1)=x(tn)+Δtdx(tn)+Δt2 d2x(tn)+Δt3 d3x(tn)+….. 1! dt 2! dt2 3! dt3
(17.3)
(17.4)
x(t )=x(t)+Δtdx(t) Δt2 d2x Δt3 d3x InitialValueProblem n+1 n 1! dt n + 2! dt2 (tn)+ 3! dt3 (tn)+…..
If we truncate Taylors series, above, and ignore terms of the order of Δt2, we will get the following approximate value for x(tn+1) (remember that tn+1 = tn + Δ),
x(tn+1) = x(tn) + Δt d x (tn) dt
x(tn + Δt) = x(tn) + Δt dx (tn) dt
(17.5) Note that the approximation above has error O(Δt2).
d =f(t, ) dt
tmin t tmax (tmin) = min
Approximate solution
Numerical/Computer Methods
Recall that ODEs are typically expressed as
dx dt
Substituting Eq. (17.6) into Eq. (17.5)
dx(tn) dt
to obtain
(17.7) where xn+1 is the approximate value of x(tn+1) and xn is the
approximate value of x(tn).The formula above is called the explicit Euler’s method and it can be used to approximate the solution at time level tn+1 when you know the (approximate) solution at time level tn.
x(tn + Δt) = x(tn) + Δt
= f(t, x) f(t,x)
(17.6)
xn+1 = xn + Δtf(tn, xn)
Smaller Δt will lead to smaller error
Predicted value of x(tn+1)
x(t)
dx
= f(tn, xn) x(tn+1)
Error
dt
Δtdx =Δtf(tn,xn) dt
x(t)
True value of x(tn+1)
x(tn)
t
tn
tn+1
xn+1 = xn + Δtf(tn, xn) Euler’s formula
Smaller Δt will lead to smaller error
(t)
dx
= f(tn, xn)
Error
dt
Predicted value of x(tn+1) x(t)
x(tn+1) x(tn)
True value of x(tn+1)
xn+1 = xn + Δtf(tn, xn) Euler’s formula
t
tn
tn+1
Smaller Δt will lead to smaller error dt
(t)
dx =f(tn,xn)
Predicted value of x(tn+1) x(tn+1)
x(tn)
x(t) True value of x(tn+1)
xn+1 = xn + Δtf(tn, xn) Euler’s formula
Error
t
tn
tn+1
x(t) xn+1 =xn+Δtf(tn,xn)=xn+Δtdx(tn,xn) dt
t
Predicted x(t) t t t
dx(t0)Δt dt
x0 = α
dx(t1)Δt dt
x2 = x1 + f (t1, x1)Δt (t2, x2)
dx(t2)Δt dt
(t1, x1)
x1 = x0 + f (t0, x0)Δt
(t3, x3)
x3 = x2 + f (t2, x2)Δt
(t0, x0)
dx(t3)Δt dt
(t4, x4) x4 = x3 + f (t3, x3)Δt
t0 = a
t1 t2 t3
t4 = b
Example L17.3:
Using Euler’s method, solve
dx dt
For 0 ≤ t ≤ 8 with x(t = 0) = 0 and a) Δt = 2
b) Δt = 1 c) Δt = 0.5 d) Δt = 0.1
=1−x
For this question f(t, x) = 1 − x
So the explicit Euler formula xn+1 = xn + Δtf(tn, xn) becomes
xn+1 = xn + Δt(1 − xn)
dx = 1 − x dt
0 ≤ t ≤ 8 x(t = 0) = 0 xn+1 = xn + Δt(1 − xn)
Lecture17M03.m
x(t) 2
Δt = 2
(t1 =0,×1 =0)
t1 = 0 t2 = 2
Δt = 2
Δt = 2
Δt = 2
(t5 =8,×5 =0)
t4 = 6 t5 = 8
(t2 = 2,×2 = 2)
(t4 = 6,×4 = 2)
(t3 =
Delta_t=2.0;
x1 = α=0
4,×3 =0)
x0 = α = 0
t3 = 4
xn+1 = xn + Δt(1 − xn)
n=1
x2 =0+2(1−0) =2
n=2
x3 =2+2(1−2) =0
n=3
x4 =0+2(1−0) =2
n=4
x5 =2+2(1−2) =0
clear all; t=0:Delta_t:8.0;
x(1)=0.0;
for n=1:length(t)-1 x(n+1)=x(n)+Delta_t*(1-x(n))
end
plot(t,x,’ko-‘,’linewidth’,2,’Markersize’,10) hold on
tplot=0:0.01:8.0; truesolution=@(t)(1-exp(-t)) plot(tplot,truesolution(tplot),’b-‘,’linewidth’,2) xlabel(‘t’);
ylabel(‘x’);
legend(‘Euler’,’True’);
clear all;
close all; Delta_t=2.0; t=0:Delta_t:8.0;
x(1)=0.0;
for n=1:length(t)-1 x(n+1)=x(n)+Delta_t*(1-x(n))
end
plot(t,x,’ko-‘,’linewidth’,2,’Markersize’,10) hold on
tplot=0:0.01:8.0; truesolution=@(t)(1-exp(-t)) plot(tplot,truesolution(tplot),’b-‘,’linewidth’,2) xlabel(‘t’);
ylabel(‘x’);
legend(‘Euler’,’True’);
Output
Δt = 2.0
Δt = 2.0 Δt = 1.0
Δt = 0.5
Δt = 0.1
End of Example L17.3
R1 = 0.5m
This example is taken from
MATLAB An Introduction by Amos Gilat.
The rate of change of the height of the water in the tank, y, is governed by
dy = − 2gyrh2 dt (2 − 0.5y)2
Initial water height is y = 2m. Find time when y = 0.1m. Given that g = 9.81m/s2.
3.0m y
r = 0.025m h
R2 = 2.0m
43
Local & Global Truncation Error
x(tn+1) = x(tn) + Δt dx (tn) 1! dt
The Taylor series estimates that the Euler method, has local truncation error O(Δt2). This means that if youhalveΔt,youcanexpecttheerrortoreducebyaquarter. ButtheTaylorseriesdoesnotgiveany information about how the error propogates, hence it does not provide any information
about the global truncation error.
However, we have to take O(Δt−1) steps to find the solution at time tN. A conservative position would to assume that the error accumulate at every time step
+Δt2 d2x(tn)+…. 2! dt2
EN = |x(tN) − xN|
global
= |x(tN) − (x0 + Δtf(t0, x0) + Δtf(t1, x1) + Δtf(t2, x2) + . . . + Δtf(tN−1, xN−1))
Hence, the global truncation error is O(Δt). In general, if a method has local truncation error of O(Δtr), then the global truncation error is O (Δtr−1).
x(t)
Δt1
t1 = b t
Δt2
t0 = a x(t)
Δt2
Error @ x = b
≈
≈ 2K((Δt1/2)2) ≈
≈ 4K((Δt1/4)2) ≈
≈ 8K((Δt1/8)2) ≈
K(Δt2) 1
Error @ x = b
K(Δt2/2) 1
t0 = a
x(t) Δt3
t0 = a
x(t) Δt4 Δt4
x1
t2 = b
t
Error @ x = b
t
Error @ x = b
Δt3
Δt3
Δt3
K(Δt2/4) 1
t1
t2
t3
t4 = b Δt4
Δt4
Δt4
Δt4
Δt4
Δt4
K(Δt2/8) 1
t0=a t1
t2
t3
t4
t5
t6
t7 t8=b
t
L17.3: Preallocation of memory
46
Preallocating Memory
Preallocating Memory
Credits:Wilson
This is a very inefficient MATLAB program!!
Lecture17M03.m
clear all;
close all;
Delta_t=2.0;
t=0:Delta_t:8.0
x(1)=0.0;
tplot=0:0.01:8.0;
truesolution=@(t)(1-exp(-t))
plot(tplot,truesolution(tplot),’b-‘,’linewidth’,2)
xlabel(‘t’);
ylabel(‘x’);
legend(‘Euler’,’True’);
for n=1:length(t)-1
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
x(n+1)=x(n)+Delta_t*(1-x(n))
end end
How to make the our Matlab program run faster: Memory Allocation
In the MATLAB program below, a new block of memory is allocated to the variable x() at every single iteration of the “for” loop
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
This is very inefficient because Matlab has to find a new block of memory for x() 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 x() before you use it in the “for” loop. You can preallocate memory and zero every element in the vector by issuing the command
x=zeros(1,100);
This command allocates memory to 100 elements in the x() vector (1 x 100) and all elements have a zero value.
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
With preallocation of memory
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
With preallocation of memory
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
With preallocation of memory
t=0:Delta_t:8.0
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
x=zeros(1,length(t));
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
phi(n+1)=phi(n)+Delta_t*(1-
phi(n))
end
With preallocation of memory
x=
0
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
x=
x=
0
0
x(2)
With preallocation of memory
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x=
x=
0
x(2)
x(3)
With preallocation of memory
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
With preallocation of memory
x=
x=
0
x(2)
x(3)
x(4)
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
x(4)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
x(4)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
x(4)
With preallocation of memory
x=
t=0:Delta_t:8.0
phi=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
x(4)
x=
x=
0
x(2)
x(3)
x(4)
x(5)
With preallocation of memory
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
x(4)
x(5)
With preallocation of memory
x=
t=0:Delta_t:8.0
x=zeros(1,length(t));
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
Without preallocation of memory
t=0:Delta_t:8.0
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
0
x(2)
x(3)
x(4)
x(5)
x=
With preallocation of memory
x=zeros(1,length(t));
0
x(02)
x(03)
x(04)
x(05)
t=0:Delta_t:8.0
for n=1:length(t)-1
x(n+1)=x(n)+Delta_t*(1-x(n))
end
x=
EXAMPLE L17.4:
Execute the MATLAB program below 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 x() 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.
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 %x=zeros(1,numel(t));
x(1)=0.0;
tic % start timer
for n=1:length(t)-1 x(n+1)=x(n)+Delta_t*(1-x(n));
end
temp=toc; %end timer
fprintf(1,’It takes %e seconds to execute the Matlab program\n’,temp);
Time it takes to execute the MATLAB program on the previous slide with and without preallocation of the array x() for different values of Δt.
WITHOUT preallocation of memory
WITH preallocation of memory
Δt=1.0e-4
3.51e-2s
2.23e-03s
Δt=1.0e-5
3.58e-1s
2.24E-02s
Δt=1.0e-6
3.81s
2.48e-1s
Δt=1.0e-7
141.1s
25.9s