Numerical Methods in Engineering (ENGR20005)
Lecture 16 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
Assignment 2
Assignment 2
Assignment 2
L16.1:
Multidimensional Boundary Value Problem
Consider the temperature distribution on a rectangular plate plate, φ(x, y), with thermal conductivity k heated by a source Q(x, y). The temperature on the plate would depend on the location of the heat
source and thermal conductivity. We would also typically know the conditions at the boundary. To make things simple, for the purpose of this lecture, we will assume that the values of φ is known
at the boundaries.
φ = φA
φ = φC
φ = φD
Q(x, y)
φ = φB
It is possible to obtain the temperature distribution on the plate by solving the 2d Poisson’s equation
22 (∂x ∂y)
k ∂2φ+∂2φ =Q(x,y)
(∂x ∂y)
I will use the notation φj,i as the approximated value of φ(xi, yj).
φj,i ≈ φ(xi, yj)
k ∂2φ+∂2φ =Q(x,y) 22
(16.1)
Please note that for 2d matrices, MATLAB stores your y data values in rows and x values in columns. Since it is common for j index to be associated with y and i index to be associated with x, I decided to use use the notation φj,i instead of φi,j.
The partial derivatives in Eq. (16.1) can be approximated using the 2nd order central differencing approximation
∂2φ (xi, yj) ≈ φj,i−1 − 2φj,i + φj,i+1 ∂x2 Δ2x
∂2φ + ≈ ∂x2
Hence Eq. (16.1) can be approximated as
φj,i−1 − 2φj,i + φj,i+1 + Δ2x
+ φj−1,i − 2φj,i + φj+1,i Δ2y
= Q(xi, yj)/k
∂2φ ∂y2
φj,i−1 − 2φj,i + φj,i+1 Δ2x
∂2φ (xi, yj) ≈ φj−1,i − 2φj,i + φj+1,i ∂y2 Δ2y
φj−1,i − 2φj,i + φj+1,i Δ2y
φj,i−1 − 2φj,i + φj,i+1 + φj−1,i − 2φj,i + φj+1,i = Q(xi, yj)/k Δ2x Δ2y
φ+φ φ+φ
⇒ j,i−1Δ j,i+1 − 2φj,i ( Δ1 + Δ1 ) + j−1,iΔ j+1,i = Q(xi, yj)/k
2x 2x2y 2y
For simplicity, I will assume that Δx = Δy = Δ. Note that this is not usually the case. I am
making this assumption to make the equations simpler and hence provide a better clarity of the solution process.
φj,i−1 + φj,i+1 + φj−1,i + φj+1,i − 4φj,i = Δ2Q(xi, yj) k
1234
(0,0)
1
2 3
4
j
φj,i−1 + φj,i+1 + φj−1,i + φj+1,i − 4φj,i =
i
v
Δ2Q(xi, yj) k
(0,0) 1 2 3 4 i
i=1,j=1:φ1,0+φ1,2+φ0,1+φ2,1−4φ1,1 = Δ2Q(x1,y1) k
i=2,j=1:φ1,1+φ1,3+φ0,2+φ2,2−4φ1,2 = Δ2Q(x2,y1) k
i=3,j=1:φ1,2+φ1,4+φ0,3+φ2,3−4φ1,3 = Δ2Q(x3,y1) k
1 2 3
4
j
i=1,j=2:φ2,0+φ2,2+φ1,1+φ3,1−4φ1,1 = Δ2Q(x1,y2) k
i=2,j=2:φ2,1+φ2,3+φ1,2+φ3,2−4φ2,2 = Δ2Q(x2,y2) k
i=3,j=2:φ2,2+φ2,4+φ1,3+φ3,3−4φ2,3 = Δ2Q(x3,y2) k
i=1,j=3:φ3,0+φ3,2+φ2,1+φ4,1−4φ3,1 = Δ2Q(x1,y3) k
i=2,j=3:φ3,1+φ3,3+φ2,2+φ4,2−4φ3,2 = Δ2Q(x2,y3) k
i=3,j=3:φ3,2+φ3,4+φ2,3+φ4,3−4φ3,3 = Δ2Q(x3,y3) k
φj,i−1 + φj,i+1 + φj−1,i + φj+1,i − 4φj,i =
Δ2Q(xi, yj) k
(16.2)
(0,0) 1 2 3 4 i
1 2 3
4
j
φ1,1 φ1,2 φ1,3 φ2,1 φ2,2 φ2,3 φ3,1 φ3,2 φ3,3
=
(16.3)
Δ2Q(x1, y1)/k − φ1,0 − φ0,1 Δ2Q(x2, y1)/k − φ0,2 Δ2Q(x3, y1)/k − φ0,3
Δ2Q(x1, y2)/k − φ2,0 Δ2Q(x2, y2)/k Δ2Q(x3, y2)/k
Δ2Q(x1, y3)/k − φ4,1 − φ3,0 Δ2Q(x2, y3)/k − φ4,2 Δ2Q(x3, y3)/k − φ4,3
−4 1 0 1 0 0 0 0 0 1 −4 1 0 1 0 0 0 0 0 1 −4 1 0 1 0 0 0
1 0 0 −4 1 0 1 0 0 0 1 0 1 −4 1 0 1 0 0 0 1 0 1 −4 1 0 1
0 0 0 1 0 0 −4 1 0 0 0 0 0 1 0 1 −4 1 0 0 0 0 0 1 0 1 −4
−41010 0 1−4101 0 01−410 1 1 00−410 0 1 0 1 −4 1
0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0
φ1,1 0 0 0 φ1,2
0 0 0 φ1,3 0 0 0 φ
1 0 0 2,1 0 1 0 φ2,2
=
Δ2Q(x1, y1)/k − φ1,0 − φ0,1
Δ2Q(x2, y1)/k − φ0,2
Δ2Q(x , y )/k − φ 31 0,3
Δ2Q(x1, y2)/k − φ2,0 Δ2Q(x2, y2)/k
Δ2Q(x3, y2)/k Δ2Q(x1, y3)/k − φ4,1 − φ3,0 Δ2Q(x2, y3)/k − φ4,2
Δ2Q(x , y )/k − φ 33 4,3
(16.3)
0 1 φ2,3 φ3,3
1 −4 1
00−41 0 φ3,1 101−41 φ 0101−4 3,2
Eφq. (16.3) is a linear algebraic equation. You should be able to solve it to obtain the values of j,i
This system of equation has lots of zeros. It is not efficient to solve lots of zeros!
You can imagine that this system quickly gets a lot bigger if we make Δ smaller and have more grid points in the domain…and you will probably have even more zeros!
It is more efficient to use Gauss-Seidel method to iteratively solve for φj,i. You can rearrange Eq. (16.2) to get
Δ2Q(x,y) φj,i=14(φj,i−1+φj,i+1+φj−1,i+φj+1,i− ki j )
(16.4)
Example L16.1: Solve Laplace’s equation ∂2φ + ∂2φ = 0 ∂x2 ∂y2
inthedomain0≤x≤4and0≤y≤3. UseΔx =Δy =Δ. Youaregiventheboundaryconditions φ(x,0) = 0oC
φ(x,3) = 30oC φ(0,y) = 20oC φ(4,y) = 10oC
y
φ(x,3) = 30oC
3
For this problem, Q(x, y) = 0. So Eq. (16.4) simplifies to
φj,i = 14 (φj,i−1 + φj,i+1 + φj−1,i + φj+1,i) φ(4,y) = 10oC
x
∂2φ + ∂2φ = 0 ∂x2 ∂y2
φ(0,y) = 20oC
(0,0)
φ(x,0) = 0oC 4
Example L16.1: Solve Laplace’s equation ∂2φ + ∂2φ = 0 ∂x2 ∂y2
inthedomain0≤x≤4and0≤y≤3. UseΔx =Δy =Δ. Youaregiventheboundaryconditions φ(x,0) = 0oC
φ(x,3) = 30oC φ(0,y) = 20oC φ(4,y) = 10oC
y
3
2 φ(0,y) = 20oC
1
φ(x,3) = 30oC
Grid for
Δx =Δy =Δ=1 φ(4,y) = 10oC
x (0,0) 1 2o 3 4
φ(x,0) = 0 C
y φ(x,3) = 30oC 301234x
NOTE that the convention in MATLAB is that increasing y in the in the direction of
increasing rows. Hence increasing y is in the downward direction.
φ(0,y) = 20oC
0
2 1
1 2
3
y(0,0) 1
φ(4,y) = 10oC
2 3 4 φ(x,0) = 0oC
x
The MATLAB commands
Delta=1.0; Lx=4.0; Ly=3.0;
x=0:Delta:Lx; y=0:Delta:Ly; nx=length(x)-1; ny=length(y)-1; [X,Y]=meshgrid(x,y);
will generate the mesh and the output below
0 1 2
3
y
01234
x
The MATLAB commands
Delta=1.0;
Lx=4.0; Ly=3.0;
x=0:Delta:Lx; y=0:Delta:Ly; nx=length(x)-1; ny=length(y)-1; [X,Y]=meshgrid(x,y);
will generate the mesh and the output below
0 1 2
3
y
NOTE that the convention in MATLAB is that increasing y in the in the direction of
increasing rows. Hence increasing y is in the downward direction.
Δx = 1 01234x
X(1,5) Y(1,5)
X(2,5) Y(2,5)
X(3,5) Y(3,5)
X(4,5) Y(4,5)
X(1,1) Y(1,1)
X(1,2) Y(1,2)
X(1,3) Y(1,3)
X(1,4) Y(1,4)
X(2,1) Y(2,1)
X(2,2) Y(2,2)
X(2,3) Y(2,3)
X(2,4) Y(2,4)
X(3,1) Y(3,1)
X(3,2) Y(3,2)
X(3,3) Y(3,3)
X(3,4) Y(3,4)
X(4,1) X(4,2) Y(4,1) Y(4,2)
X(4,3) X(4,4) Y(4,3) Y(4,4)
Δy = 1
The MATLAB command
phi=zeros(size(X));
will define the matrix φj,i which is the approximation of φ(xi, yj) and set its values to 0.
01234
0 1 2
3
y
x
φ(1,1) =0
φ(2,1) =0
φ(1,2) =0
φ(3,2) =0
φ(1,3) =0
φ(3,3) =0
φ(1,4) =0
φ(3,4) =0
φ(1,5) =0
φ(2,2) =0
φ(2,3) =0
φ(2,4) =0
φ(2,5) =0
φ(3,1) =0
φ(3,5) =0
φ(4,1) =0
φ(4,2) =0
φ(4,3) =0
φ(4,4) =0
φ(4,5) =0
The MATLAB commands
Ttop=30.0; Tbottom=0.0; Tleft=20; Tright=10;
Sets up the boundary conditions
0 1 2
3
y
φ(x,0) = 0oC 01234
x
phi(:,1)=Tleft;
phi(:,end)=Tright;phi(1,:)=Tbottom;phi(end,:)=Ttop;
φ(1,1) =0
φ(2,1) = 20
φ(1,2) =0
φ(2,2) =0
φ(1,3) =0
φ(2,3) =0
φ(1,4) =0
φ(2,4) =0
φ(1,5) =0
φ(2,5) = 10
φ(3,1) = 20
φ(4,1) = 30
φ(3,2) =0
φ(4,2) = 30
φ(3,3) =0
φ(4,3) = 30
φ(3,4) =0
φ(3,5) = 10
φ(4,4) = 30
φ(4,5) = 30
φ(x,3) = 30oC
φ(0,y) = 20oC
φ(4,y) = 10oC
The MATLAB commands
%%
avgerror=1;
while avgerror>1.0e-6 for j=2:ny
for i=2:nx phi(j,i)=(phi(j-1,i)+phi(j+1,i)+phi(j,i-1)+phi(j,i+1))/4;
end end
avgerror=0.0;
for j=2:ny for i=2:nx
avgerror=avgerror+abs(4*phi(j,i)-(phi(j-1,i)+phi(j+1,i)+phi(j,i-1)+phi(j,i+1)));
end end
avgerror/(nx*ny)
end
0 until avgerror is smaller than 10−6.
Using Gauss-Seidel (see L6.5) to solve
a system of equation similar to Eq. (16.4) using
φj,i = 14 (φj,i−1 + φj,i+1 + φj−1,i + φj+1,i)
φ(x,0) = 0oC 01234
x
iteratively calculates the value of φj,i
1
φ(1,1) =0
φ(2,1) = 20
φ(1,2) =0
φ(2,2) = 12.75
φ(1,3) =0
φ(2,3) = 10.43
φ(1,4) =0
φ(2,4) = 9.42
φ(1,5) =0
φ(2,5) = 10
φ(3,1) = 20
φ(4,1) = 30
φ(3,2) = 20.58
φ(4,2) = 30
φ(3,3) = 19.57
φ(4,3) = 30
φ(3,4) = 17.25
φ(3,5) = 10
φ(4,4) = 30
φ(4,5) = 30
2
3
y
φ(x,3) = 30oC
φ(0,y) = 20oC
φ(4,y) = 10oC
The MATLAB commands
[C,h]=contour(X,Y,phi,10,’linewidth’,5) clabel(C,h)
h.LevelList=round(h.LevelList,2) xlabel(‘x’,’Fontsize’,14);ylabel(‘y’,’Fontsize’,14);
produces the contour plot of φj,i ≈ φ(xi, yj)
φ(x,0) = 0oC 01234x
0 1 2
3
y
φ(4,1) = 30
φ(4,2) = 30
φ(4,3) = 30
φ(4,4) = 30
φ(1,5) =0
φ(2,5) = 10
φ(3,5) = 10
φ(4,5) = 30
φ(1,1) =0
φ(1,2) =0
φ(1,3) =0
φ(1,4) =0
φ(2,1) = 20
φ(2,2) = 12.75
φ(2,3) = 10.43
φ(2,4) = 9.42
φ(3,1) = 20
φ(3,2) = 20.58
φ(3,3) = 19.57
φ(3,4) = 17.25
φ(x,3) = 30oC
φ(0,y) = 20oC
φ(4,y) = 10oC
Running the code again with smaller values of Δ will produce the plot below
φ(x,3) = 30oC
φ(x,0) = 0oC
φ(0,y) = 20oC
φ(4,y) = 10oC
clear all
Delta=0.1; Lx=4.0; Ly=3.0;
Ttop=30.0; Tbottom=0.0; Tleft=20; Tright=10;
%% x=0:Delta:Lx; y=0:Delta:Ly; nx=length(x)-1; ny=length(y)-1;
[X,Y]=meshgrid(x,y); phi=zeros(size(X));
% Setting the Dirichlet Boundary Conditions
Ttop=30.0;
Tbottom=0.0;
Tleft=20;
Tright=10; phi(:,1)=Tleft;phi(:,end)=Tright;phi(1,:)=Tbottom;phi(end,:)=Ttop; %%
avgerror=1;
while avgerror>1.0e-6 for j=2:ny
for i=2:nx phi(j,i)=(phi(j-1,i)+phi(j+1,i)+phi(j,i-1)+phi(j,i+1))/4;
end end
avgerror=0.0;
for j=2:ny for i=2:nx
avgerror=avgerror+abs(4*phi(j,i)-(phi(j-1,i)+phi(j+1,i)+phi(j,i-1)+phi(j,i+1)));
end end
avgerror/(nx*ny)
end
[C,h]=contour(X,Y,phi,10,’linewidth’,5) clabel(C,h)
h.LevelList=round(h.LevelList,2) xlabel(‘x’,’Fontsize’,14);ylabel(‘y’,’Fontsize’,14);
φ(x,3) = 30oC
φ(x,0) = 0oC
Lecture16M01.m
φ(0,y) = 20oC
φ(4,y) = 10oC
Example L16.1: Solve Laplace’s equation ∂2φ + ∂2φ = 0 ∂x2 ∂y2
inthedomain0≤x≤4and0≤y≤3. UseΔx =Δy =Δ. Youaregiventheboundaryconditions
Mixed BC
∂φ (x,0) = 0 ∂y o
φ(x,3) = 30 C φ(0,y) = 20oC φ(4,y) = 10oC
y
φ(x,3) = 30oC
3
For this problem, Q(x, y) = 0. So Eq. (16.4) simplifies to
φj,i = 14 (φj,i−1 + φj,i+1 + φj−1,i + φj+1,i) φ(4,y) = 10oC
∂2φ + ∂2φ = 0 ∂x2 ∂y2
φ(0,y) = 20oC
(0,0)
4
x
∂φ (x,0) = 0 ∂y
L16.2:
Spectral differentiation for
intervals other than −1 ≤ x ≤ 1
25
In lecture 14, we introduced the concept of the derivative matrix. We can use this matrix to very accurately evaluate the derivative of any
function, f(x).
Derivative matrix
L0′(x0) L1′(x0) L2′(x0)
L0′(x1) L1′(x1) L2′(x1)
. . . . . . . . .
Ln′(x0) Ln′(x1) Ln′(x2)
L0′(x2) L1′(x2) L2′(x2)
. . …..
. . ….. L0′(xn) L1′(xn) L2′(xn) . . . Ln′(xn)
fn′(x0) fn′(x1) =
fn′(x2) .
f0
f1
f2 .
fn
accuracy for the fi′ if the derivative matrix is evaluated at x = xi where xi are the
roots of the Legendre or Lobatto polynomials (GLe and GLo points). These
functions are usually defined between -1 and 1 and hence the derivative matrix defined
above is only really valid in the domain −1 ≤ x ≤ 1. 26
fn′(xn)
Li(x) are your Legendre basis functions and I showed that you can get very good
If we would like to differentiate a function, f(t), for a ≤ t ≤ b, we can transform this Interval to −1 ≤ x ≤ 1 using just a simple linear transformation
t = px + q
We need to find p and q that relates x and t. We would like
x=−1 when x=+1
27
t=a when t=b
x=−1 when t=a x=+1 when t=b
a=−p+q b=p+q
We can solve the above two equations to get
t = px + q Hence
p=b−a 2
q=b+a 2
Therefore
t=b−ax+b+a 22
x= 2 t−b+a (b−a) b−a
dx= 2
dt (b−a) 28
x= 2 t−b+a (b−a) b−a
df = df dt dx
=( 2 )df b−a dx
dx= 2
dt (b−a)
dx dt
df dt
So to calculate df/dt we just need to calculate df/dx and multiply by 2/(b − a)
It can also be shown that 2 ( ) df=d df
dt2 dt
d2f=d( 2 df)
dt
dt2 dt b − a dx
d2f=( 2 )2d2f
2 2 dt2 b−a dx2 2 2
So to calculate d f/dt we just need to calculate d f/dx and multiply by
4/(b − a)2.
29
Example L16.2: Write a MATLAB program that implements Eq. (14.1) to approximate the derivative of
f(t) = et sin(t)
In the domain t ∈ [−3,2]. Use the GLo points. Compare with the exact answer
df = et(sin(5t) + 5 cos(5t)) dt
30
Lecture16M02.m
tpoints=((b-a)/2)*xpoints+(a+b)/2;
fpoints=exp(tpoints).*sin(5*tpoints);
t=b−ax+b+a 22
f(t) = et sin(t)
derivfpoints=(2/(b-a))*D*fpoints;
df = et(sin(5t) + 5 cos(5t)) dt
dfreal=exp(tint).*(sin(5*tint)+5*cos(5*tint));
df=( 2 )df dt b−a dx
clear all a=-3; b=2;
%xpoints=linspace(a,b,11)’; %[xpoints,w]=lgwt(21,-1,1); [xpoints,w,P]=lglnodes(51);
n=length(xpoints)-1; %order of polynomial
xint=-1:0.01:1;
tint=((b-a)/2)*xint+(a+b)/2; fint=LagrangePolynomial(tpoints,fpoints,n,tint); freal=exp(tint).*sin(5*tint);
figure(1)
hold off plot(tpoints,fpoints,’ko’,’MarkerSize’,20,’MarkerFaceColor’,’r’) hold on
plot(tint,fint,’b-‘,’Linewidth’,4)
plot(tint,freal,’k-‘,’Linewidth’,4)
xlabel(‘t’);ylabel(‘f(t)’)
D=DerivMatrix(xpoints,n);
figure(2)
hold off
plot(tint,dfreal,’k-‘,’Linewidth’,4)
hold on plot(tpoints,derivfpoints,’ko’,’MarkerSize’,20,’MarkerFaceColor’,’b’) xlabel(‘t’);ylabel(‘df(t)/dt’)
31
In lecture 14, we introduced the concept of the derivative matrix. We can use this matrix to very accurately evaluate the derivative of any function
We can also view the second derivative as a matrix operation.
d2f(xi)= d (df)(xi) dx2 dx dx
≈ [D]([D]{f}) = [D][D]{f} = [D]2{f}
So the second derivative matrix can be approximated as the square of the first derivative matrix. Only for spectral differentiation!
32
L16.3:
Solving Boundary Value Problem Using Spectral Methods
33
In lecture 14, we introduced the concept of the derivative matrix where elements of the matrix are the derivative of the Lagrange basis functions evaluated at equally spaced, GLo and GLe points. We can use this matrix to very accurately evaluate the derivative of any
function, f(x). In this section, we are going to see how these matrices can be used to accurately solved the boundary value problem.
It is probably easiest to illustrate the method by doing an example.
Consider a well posed one-dimensional two point boundary value problem of the form
y′′ = e4x with dirichlet boundary conditions y(−1) = α and y(1) = β.
(16.5) For the moment, we will restrict our analysis and only consider the domain −1 ≤ x ≤ 1.
Discretising the second derivative in Eq. (16.5) with the second derivative matrix gives the linear system
[D2]{y} = {f}
Since the boundary conditions are inhomogeneous, we’ll look at the structure of the system to get an idea of how these interact with the interior nodes. To simplify the analysis, and to get a clearer idea on how things work, we will only consider the case where there are 5 intervals in −1 ≤ x ≤ 1.
34
y(−1) = α
y0 = α y1 y2
x0=−1 x1 x2
y(1) = β y5 = β
y3 y4
x3 x4 x5=1
y′′ = e4x =
y′′ 0
y′′ 1
D2 D2 D2 D2 D2 D2 00 01 02 03 04 05
y D2 D2 D2 D2 D2 D2 0
y
D D D 20 21 22
D2 D2 D2 30 31 32
D2 D2 D2 40 41 42
D2 D2 D2 50 51 52
D D D y 23 24 25 2
y′′ 3
y′′ 4
y′′ 5
y 4
′′
10 11 12 13 14 15 y 2 2 2 2 2 2 1
2
=
D2 D2 D2 y 33 34 35 3
D2 D2 D2 43 44 45
D2 D2 D2 53 54 55
y 5
35
e4x0 e4x1 e4x2 e4x3 e4x4 e4x5
y(−1) = α y0 = α
y(1) = β y5 = β
y1y2 y3y4
D2 D2 D2 D2 D2 D2 10 11 12 13 14 15
D2 D2 D2 D2 D2 D2 20 21 22 23 24 25
D2 D2 D2 D2 D2 D2 30 31 32 33 34 35
D2 D2 D2 D2 D2 D2 40 41 42 43 44 45
e4x1 e4x2 e4x3 e4x4
100000
D2 D2 D2 D2 D2 D2 10 11 12 13 14 15
D2 D2 D2 D2 D2 D2 20 21 22 23 24 25
D2 D2 D2 D2 D2 D2 30 31 32 33 34 35
D2 D2 D2 D2 D2 D2 40 41 42 43 44 45
000001
e4x1 e4x2 e4x3 e4x4
y′′ = e4x x0: D2 D2 D2 D2 D2 D2
x3 =
e4x5 α
=
β
x4 x5=1
x0=−1 x1 x2
x1: 00 01 02 03 04 05
x2 : x3 :
y0 y1 y2 y3 y4 y5
y0 y1 y2 y3 y4 y5
e4x0
x4 :
x5: D2 D2 D2 D2 D2 D2
50 51 52 53 54 55
36
y0 = α
y1y2 y3y4
y5 = β
D2 D2 10 11
D2 D2 20 21
D2 D2 30 31
D2 D2 40 41
D2 D2 12 13
D2 D2 22 23
D2 D2 32 33
D2 D2 42 43
D2 D2 14 15
D2 D2 24 25
D2 D2 34 35
D2 D2 44 45
x0=−1 x1
x2
D2 D2 D2 D2 D2 D2
x3 x4 x5=1 e4x0
e4x1
= e4x2 e4x3 e4x4 e4x5
00 01 02 03 04 05
y0 y1 y2 y3 y4 y5
D2 D2 D2 D2 D2 D2 50 51 52 53 54 55
37
y0 = α
y1y2 y3y4
y5 = β
D2 D2 D2 D2 D2 D2 10 11 12 13 14 15
D2 D2 D2 D2 D2 D2 20 21 22 23 24 25
D2 D2 D2 D2 D2 D2 30 31 32 33 34 35
D2 D2 D2 D2 D2 D2 40 41 42 43 44 45
x0=−1×1 x2 x3 x4
D2 D2 D2 D2 D2 D2 00 01 02 03 04 05
x5 = 1
e4x0
y0
y1
y2 = e4x2
D2 D2 D2 D2 D2 D2 50 51 52 53 54 55
D2 D2 D2 D2 D2 D2 10 11 12 13 14 15
D2 D2 D2 D2 D2 D2 20 21 22 23 24 25
D2 D2 D2 D2 D2 D2 30 31 32 33 34 35
D2 D2 D2 D2 D2 D2 40 41 42 43 44 45
y3 y4 y5
e4x3 e4x4 e4x5
y1 e4x1
e4x1
y0
38
y 4x 2 =e2
y3 e4x3
y4 e4x y 4
5
y0 = α
y1y2 y3y4
y5 = β
D2 10
D2 20
D2 30
D2 40
D2 15
D2 25
D2 35
D2 45
D2 10
D2 20
D2 30
D2 40
D2 15
D2 25
D2 35
D2 45
x0=−1 x1 x2 x3 y
x4 x5=1
D2 D2 D2 D2 y0 e4x 11 12 13 14 1 1
D2 D2 D2 D2 y 4x 21 22 23 24 2 =e2
D2 D2 D2 D2 y3 e4x3 31 32 33 34 y4 e4x
D2 D2 D2 D2 y 4 41 42 43 44 5
D2 D2 D2 D2 y 11 12 13 14 1
e4x 1
4x 2+y =e2
D2 D2 D2 D2 y
y +21222324 0 D2 D2 D2 D2
y3 5 e4x3
31 32 33 34
y4
e4x 4
D2 D2 D2 D2 41 42 43 44
39
y0 = α
y1y2 y3y4
y5 = β
D2 10
D2 20
D2 30
D2 40
D2 15
D2 25
D2 35
D2 45
D2 10
D2 20
D2 30
D2 40
D2 15
D2 25
D2 35
D2 45
x0=−1 x1 x2 x3 x4 x5=1 D2 D2 D2 D2 y e4x
11 12 13 14 1 1 D2 D2 D2 D2 y 4x
Rearranging gives
y +21222324 2+y =e2 0 D2 D2 D2 D2 y3 5 e4x3
D2 D2 D2 D2 y 11 12 13 14 1
e4x 1
31 32 33 34 y4 D2 D2 D2 D2
e4x 4
41 42 43 44
D2 D2 D2 D2 y
D2 D2 D2 D2 y3 e4x3
4x
21222324 2=e2−α −β
31 32 33 34 y4
D2 D2 D2 D2 4
41 42 43 44
e4x
Can solve to find y1, y2, y3 and y4. 40
y0 = α
y1 y2 yn−2 yn−1
yn = β
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
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 n−1,0
41
Lecture16M03.m
y′′ = e4x
y(−1) = α
y(1) = β
y0 y1 y2
y⋮ n−1
yn
y0 y1 y2
y⋮ n−1
yn
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
D=DerivMatrix(xpoints,n);
D2 = D*D;
Amatrix=D2(2:end-1,2:end-1);
Cvector=r(2:end-1)-alpha*D2(2:end-1,1)-beta*D2(2:end-1,end);
e4x D2
1 10 D2
4×2 D2 1n e −y 20 −y D2
⋮ 0 ⋮ n 2n e4x ⋮ D2
n−1 D2 n−1,n n−1,0
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
clear all
Npoints=50;
alpha = 1; Beta = 0;
%BC1 y(-1)=1 %BC2 y(1)=0
y0′
y1′
%xpoints=linspace(-1,1,Npoints)’;
[xpoints,w,P]=lglnodes(Npoints); xpoints=-xpoints;
= y n′ − 1
yn′ y = @(x) (exp(4*x)-sinh(4)*x-cosh(4))./16+1-0.5*(1+x);
y′′ 0
y′′ 1
y2′ ⋮
yexact = y(xpoints);
n=length(xpoints)-1; %order of polynomial
r = exp(4*xpoints);
sol = Amatrix\Cvector; yapprox = [alpha;sol;beta];
hold off plot(xpoints,yapprox,’bo’) hold on plot(xpoints,yexact,’k-‘)
function D=DerivMatrix(x,n) D=zeros(n+1,n+1);
y′′ ⋮2
= y′′
n−1 y′′
n
y1 y2
y⋮ n−1
=
42
[xpoints,w,P]=lglnodes(Npoints); xpoints=-xpoints;
y = @(x) (exp(4*x)-sinh(4)*x-cosh(4))./16+1-0.5*(1+x); yexact = y(xpoints);
n=length(xpoints)-1; %order of polynomial r = exp(4*xpoints);
sol = Amatrix\Cvector; yapprox = [alpha;sol;beta];
hold off plot(xpoints,yapprox,’bo’) hold on plot(xpoints,yexact,’k-‘)
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;
43
D=DerivMatrix(xpoints,n);
D2 = D*D;
Amatrix=D2(2:end-1,2:end-1);
Cvector=r(2:end-1)-alpha*D2(2:end-1,1)-beta*D2(2:end-1,end);
Consider now a well posed one-dimensional two point boundary value problem of the form
y′′ = e4x
Discretising the second derivative in Eq. (16.5) with the second derivative matrix gives the linear system
(16.5)
with dirichlet boundary conditions y(−1) = 1 but with Neumann boundary condition y′(1) = 0.
[D2]{y} = {f}
The difference between this equation system of equation and the previous one is the boundary condition. We’ll look at the structure of the system to get an idea of how these interact with the interior nodes. Like the last problem we considered, we will consider the problem first for only 6 grid points.
44
y(−1) = 1
y0 = α y1 y2
x0=−1 x1 x2
y3 y4 x3 x4
= e4x =
y′(1) = 0 y′ = β
5
e4x0 e4x1 e4x2 e4x3 e4x4 e4x5
x5=1
y′′
y′′ 0
y′′ 1
D2 D2 D2 D2 D2 D2 00 01 02 03 04 05
y D2 D2 D2 D2 D2 D2 0
y
D D D 20 21 22
D2 D2 D2 30 31 32
D2 D2 D2 40 41 42
D2 D2 D2 50 51 52
D D D y 23 24 25 2
y′′ 3
y′′ 4
y′′ 5
y 4
′′
10 11 12 13 14 15 y 2 2 2 2 2 2 1
2
=
D2 D2 D2 y 33 34 35 3
D2 D2 D2 43 44 45
D2 D2 D2 53 54 55
y 5
45
y1y2 y3y4
y′ = β 5
y0 = α
D2 10
D2 20
D2 30
D2 40
D2 50
D2 D2 11 12
D2 D2 21 22
D2 D2 31 32
D2 D2 41 42
D2 D2 51 52
D2 D2 D2 13 14 15
D2 D2 D2 23 24 25
D2 D2 D2 33 34 35
D2 D2 D2 43 44 45
D2 D2 D2 53 54 55
e4x 1
e4x2 e4x3
e4x 4
e4x 5
D2 D2 D2 D2 D2 D2 10 11 12 13 14 15
D2 D2 D2 D2 D2 D2 20 21 22 23 24 25
D2 D2 D2 D2 D2 D2 30 31 32 33 34 35
D2 D2 D2 D2 D2 D2 40 41 42 43 44 45
D2 D2 D2 D2 D2 D2 50 51 52 53 54 55
e4x1 e4x2 e4x3 e4x4 e4x5
x0=−1 x1
Ignore the first row since we know the value of y0
x0 : x:
00 01 02 03 04 05
y1
x2
D2 D2 D2 D2 D2 D2
x3 =
y0 y1 y2 y3 y4 y5
x4 x5=1
e4x0
y2 y3 y4 y5
x1 : 2
y0
x3 : x4 : x5 :
46
=
y1y2 y3y4
y′ = β 5
y0 = α
D2 10
D2 20
D2 30
D2 40
D2 50
D2 10
D2 20
D2 30
D2 40
D2 50
x0=−1 x1
x2 x3 x4 x5=1
e4x1
D2 D2 D2 D2 D2 y 1112131415 0
D2 D2 D2 D2 D2 y1 2122232425 y2
e4x e4x3
e4x 4
4x e5
D2 D2 D2 D2 D2 y= 3132333435 3
D2 D2 D2 D2 D2 y 4142434445 y4
D2 D2 D2 D2 D2 5 51 52 53 54 55
D2 D2 D2 D2 D2 1112131415 y1
D2 D2 D2 D2 D2 2122232425 y2
D2 D2 D2 D2 D2 y3 3132333435 y
D2 D2 D2 D2 D2 4 4142434445 y5
D2 D2 D2 D2 D2 51 52 53 54 55
e4x1 4x
2
y+ 0
e2 = e4x3
47
e4x 4
4x e5
y0 = α
y1y2 y3y4
y′ = β 5
D2 10
D2 20
D2 30
D2 40
D2 50
y5′
D50 D51 D52 D53 D54 D55
y0 y1 y2 y3 y4 y5
D2 10
D2 20
D2 30
D2 40
D2 50
x0=−1 x1
D2 D2 D2 D2
x2 15 y1
x3 x4 D2 D2 D2 D2 D2
x5=1 −y
11 12 13 14
D2 D2
e4x1 4x
1112131415 y1 ⇒D2 D2 D2 D2 D2
e4x1 4x
y+ 0
D2 D2 D2 D2 21 22 23 24
D2 D2 D2 D2 31 32 33 34
25 y2 D2 y3
=
e2
2122232425 y2 D2 D2 D2 D2 D2 y3
e2 = e4x3
D2 D2 D2 D2 41 42 43 44
D2 4
35 y
e e
3132333435 y D421 D422 D423 D424 D425 y4
e4x 4
0
45 y5 D2
22222 5 D51 D52 D53 D54 D55
4×3 e4x4
4×5
We also know that the derivative of y is given by
y0′ D00 D01 D02 D03 D04 D05
y1′ D10 D11 D12 D13 D14 D15
y4′ D40 D41 D42 D43 D44 D45 y1
D2 D2 D2 D2
51 52 53 54 55
e4x5
y5′ = D50y0 + D51y1 + D52y2 + D53y3 + D54y4 + D55y5
y2′ D20 D21 D22 D23 D24 D25
= Dy+Dy+Dy+Dy+Dy+Dy=β
y′ D D D D D D y 3 30 31 32 33 34 35
Since we are given that y′ = β, so 5
50 0 51 1 52 2 53 3 54 4 55 5
48
D50y0 + [D51 D52 D53 D54 D55]
2
y3 = β
y4 5
y0 = α
y1y2 y3y4
y′ = β 5
D2 D2 D2 D2 D2 11 12 13 14 15
D2 D2 D2 D2 D2 21 22 23 24 25
D2 D2 D2 D2 D2 31 32 33 34 35
D2 D2 D2 D2 D2 41 42 43 44 45
e4x1 e4x2 e4x3 e4x4
D2 10
D2 20
D2 30
D2 40
D2 D2 D2 D2 D2 11 12 13 14 15
D2 D2 D2 D2 D2 21 22 23 24 25
D2 D2 D2 D2 D2 31 32 33 34 35
D2 D2 D2 D2 D2 41 42 43 44 45
D51 D52 D53 D54 D55
y1 Dy+D D D D D y3
50 0 [ 51 52 53 54 55] y2 y4
=β
y5
e4x1 e4x2 e4x3 e4x4
y1 y2
[D51 D52 D53 D54 D55] y3 = β − y0D50 y4
y5
x0=−1 x1 D2 D2 D2 D2 D2
x2 x3 D2
D2 20
x4 x5=1 = −y0
β
y1 y2
10
51 52 53 54 55
e5 D2 50
y3= −yD2 y 030
4 D2 y5 4x 40
y1 y2 y3 y4 y5
D50
Rearranging gives
49
y0 = α
y1 y2 yn−2 yn−1
y′ = β n
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 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
50
Lecture16M04.m
y′′ = e4x
clear all Npoints=30;
alpha = 1; beta = 0;
y(−1) = α
y′(1) = β
y0 y1 y2
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
D=DerivMatrix(xpoints,n);
D2 = D*D;
%BC1 y(-1)=alpha %BC2 y'(1)=beta
y0′
y1′
y2′ ⋮
yn′ −1 yn′
=
y′′
%xpoints=linspace(-1,1,Npoints)’;
[xpoints,w,P]=lglnodes(Npoints); xpoints=-xpoints;
y ⋮ n−1
y = @(x) -(1-16*exp(4)+4*exp(8)-exp(4+4*x)+4*exp(8)*x)./ (16*exp(4));
Amatrix=[D2(2:end-1,2:end);D(end,2:end)];
Cvector=[r(2:end-1);beta]-alpha*[D2(2:end-1,1);D(end,1)];
y′′ n−1
yn−1 yn
y′′ n
yn
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
yexact = y(xpoints);
n=length(xpoints)-1; %order of polynomial
r = exp(4*xpoints);
sol = Amatrix\Cvector; yapprox = [alpha;sol];
hold off plot(xpoints,yapprox,’bo’) hold on plot(xpoints,yexact,’k-‘) xlabel(‘x’);ylabel(‘y’)
function D=DerivMatrix(x,n) D=zeros(n+1,n+1);
0 y0 y′′
1 y1 y′′ y2
⋮2= ⋮
D2 D2 D2 … D2 11 12 13 1,n
D2 D2 D2 … D2 21 22 23 2,n
D2 D2 D2 … ⋮ 31 32 33
D2 D2 D2 … D2
41 42 43 n−1,n
Dn,1 Dn,2 Dn,3 … Dn,n
e4x D2 1 10
y1
y2
y⋮ =
n−1 yn
51
4×2 D2 e 20
⋮−α⋮ e4xn−1 D2
β n−1,0 Dn0
for i=1:n+1
D=DerivMatrix(xpoints,n);
D2 = D*D;
Amatrix=[D2(2:end-1,2:end);D(end,2:end)];
Cvector=[r(2:end-1);beta]-alpha*[D2(2:end-1,1);D(end,1)];
n=length(xpoints)-1; %order of polynomial r = exp(4*xpoints);
sol = Amatrix\Cvector; yapprox = [alpha;sol];
hold off plot(xpoints,yapprox,’bo’) hold on plot(xpoints,yexact,’k-‘) xlabel(‘x’);ylabel(‘y’)
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