Numerical Methods in Engineering (ENGR20005)
Lecture 15 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L15.1:
Boundary Value Problem:
Finite Difference Approach (Introduction)
2
Lecture Notes (Chap. 6, pg. 105)
Boundary Value Problem
In many engineering problems, you will have to solve ordinary differential equations that look something like
where
We are also usually given the conditions at the boundaries
d2y + p(x) dy + q(x)y = r(x) dx2 dx
(15.1)
a≤x≤b y(a) = α
y(b) = β
Solution to Eq. (15.1)?
Solution to Eq. (15.1)?
a
Values of y given at the boundaries: Dirichlet boundary conditions.
y
α β
Solution to Eq. (15.1)?
x b
3
Boundary Value Problem
In many engineering problems, you will have to solve ordinary differential equations that look something like
d2y + p(x) dy + q(x)y = r(x) dx2 dx
(15.1)
Value of derivative
y given at the x = a: Neumann boundary condition. Value of y given at the x = b : Dirichlet boundary condition
Solution to Eq. (15.1)?
Solution to Eq. (15.1)?
where
We are also usually given the conditions at the boundaries
y
a≤x≤b y′(a) = α
y(b) = β Solution to Eq. (15.1)?
a
α β
x b
α α
4
Simple Beam Problem
A simply supported beam with a modulus of elasticity E, moment of inertia I, a uniform load w and end tension T is shown below.
w TT
The deflection at each point x in the beam can be calculated by solving the following differential equation:
q(x) = − a
y
x
Lecture recap
y ′′ − a y = b x ( L − x ) where L is the length of the beam and,
+ p(x) dy + q(x)y = r(x)
a=T,b=w
0
dx r(x) = bx(L − x)
EI
2EI
5
d2y dx2
=
The idea!
1. Approximate the derivatives by Finite Difference Approximation. I will use 2nd order central schemes. But you can use any finite difference approximation taught in Lecture 13.
2. Implement boundary conditions. y(a) = α Dirichlet
y′(a) = α Neumann at x = a and y(b)=β Dirichletatx=b
[A]{X} = {C} {f} = {0}
y(b)=β
3. Solve the resulting set of algebraic equations.
d2y + p(x) dy + q(x)y = r(x) dx2 dx
d2y ≈ yi−1 −2yi +yi+1 dy ≈ yi+1 −yi−1 dx2Δ2 dx2Δ
6
L15.2:
Boundary Value Problem: Finite Difference Approach (Dirichlet Boundary Conditions)
7
Lecture Notes (Chap. 6, pg. 105)
First consider the case where you “discretise” your domain a ≤ x ≤ b into 5 intervals.
(15.1)
Values of y given at the boundaries: Dirichlet boundary conditions.
d2y + p(x) dy + q(x)y = r(x) dx2 dx
y(a)=y0 =α
y(b)=y5 =β
y0
y1 y2 y3 y4
y5
ΔΔΔΔΔ
x0 x1 x2 x3 x4 x5
y0 = α
So you need to solve for the unknowns y1, y2, y3 and y4. 8
y′′ + p(x )y′ + 111
q(x1)y1 = r(x1)
y′′ + p(x )y′ + 222
q(x2)y2 = r(x2)
y′′ + p(x )y′ + 333
q(x3)y3 = r(x3)
y′′ + p(x )y′ + 444
q(x4)y4 = r(x4)
y5 = β
y0 y1 y2 y3 y4 y5 ΔΔΔΔΔ
x0 x1 x2 x3 x4 x5
y0 = α
Using
Central difference approximation
q(x2)y2 = r(x2) q(x3)y3 = r(x3) ++=
Using Central difference approximation
y5 = β
y′′ + p(x )y′+ y′′ + p(x )y′+ y′′ + p(x )y′+ 111222333
y′′ + p(x )y′ + 444
q(x4)y4 = r(x4)
q(x1)y1 = r(x1)
y′′ 1
y′′ 2
y′′ 3
y′′ 4
p(x )y′ 1 1
p(x )y′ 2 2
p(x )y′ 3 3
p(x )y′ 4 4
q(x )y 1 1
q(x )y 2 2
q(x )y 3 3
q(x )y 4 4
r(x ) 1
r(x ) 2
r(x ) 3
r(x ) 4
y1′′ = y2 − 2y1 + y0 Δ2
y2′′ = y3 − 2y2 + y1 Δ2
y3′′ = y4 − 2y3 + y2 Δ2
y4′′ = y5 − 2y4 + y3 Δ2
q(x1)y1 q(x2)y2 q(x3)y3 q(x4)y4
r(x1) r(x2) r(x3) r(x4)
9
p(x1)y1′ = p(x1) y2 − y0 2Δ
p(x2)y2′ = p(x2) y3 − y1 2Δ
p(x3)y3′ = p(x3) y4 − y2 2Δ
p(x4)y4′ = p(x4) y5 − y3 2Δ
y0 y1 y2 y3 y4 y5 ΔΔΔΔΔ
x0 x1 x2 x3 ′′x4′ x5
y =α 0
y′′+p(x)y′+ y′′+p(x)y′+ y′′+p(x)y′+ y4+p(x4)y4+ y =β 1 11 2 22 3 33 q(x)y=r 5
q(x1)y1 = r1 q(x2)y2 = r2 q(x3)y3 = r3 4 4 4 ++=
r(x1) r(x2) r(x3) r(x4)
q(x1)y1 q(x2)y2 q(x3)y3 q(x4)y4
r(x1) r(x2) r(x3) r(x4)
p(x1)y1′ p(x2)y2′ p(x3)y3′ p(x4)y4′
p ( x 1 ) y 1′ = p ( x 1 ) y 2 − y 0 2Δ
p ( x 2 ) y 2′ = p ( x 2 ) y 3 − y 1 2Δ
p ( x 3 ) y 3′ = p ( x 3 ) y 4 − y 2 2Δ
p ( x 4 ) y 4′ = p ( x 4 ) y 5 − y 3 2Δ
− 12 p ( x 1 ) 0 12 p ( x 1 ) 0 0 0 y0 0 −1p(x) 0 1p(x) 0 0 y1
12222 y2
Δ 0 0 −1p(x3) 0 1p(x3) 0 y3 2 1 2 1 y4 0 0 0 −2 p(x4) 0 2 p(x4) y5
y′′ 1
y′′ 2
y′′ 3
y′′ 4
y1′′ = y2 − 2y1 + y0 Δ2
y2′′ = y3 − 2y2 + y1 Δ2
y3′′ = y4 − 2y3 + y2 Δ2
y4′′ = y5 − 2y4 + y3 Δ2
y0 1 −2 1 0 0 0 y1
1 0 1 −2 1 0 0 y2 Δ2 0 0 1 −2 1 0 y3
0 0 0 1 −2 1 y4 y5
++=
10
q(x1)y1 q(x2)y2 q(x3)y3 q(x4)y4
0q(x) 0 0 0 0 y0 1 y1
0 0 q(x2) 0 0 0 y2
r(x1) r(x2) r (x3) r(x4)
0 0 0 q(x3) 0 0 y3
0 0 0 0 q(x4)0 y4 y5
y0 1 −2 1 0 0 0 y1
−1p(x) 0 1p(x) 0 0 0 21 21
+ 1 0 − 12 p ( x 2 ) 0 12 p ( x 2 ) 0 0
y0 y1 y 2 y3 y4 y5
y0 0q(x1) 0 0 0 0 y1
+ 0 0 q ( x 2 ) 0 0 0 y 2
r(x ) 1
= r ( x 2 ) r (x3) r(x4)
1 0 1 − 2 1 0 0 y 2 Δ2 0 0 1 −2 1 0 y3
Δ 0 0
0 −1p(x) 0 1p(x) 0 23 23
0 0 − 12 p ( x 4 ) 0 12 p ( x 4 )
0 0 0 q(x3) 0 0 y3
0 0 0 1 −2 1 y4 y5
0 0 0 0 q(x)0 y4 4 y5
y0 = α
y5 = β
100000
1 −p(x1) −2 +q(x1) 1 +p(x1) Δ2 2Δ Δ2 Δ2 2Δ
0 0 0
1 +p(x2) 0 0 Δ2 2Δ
0 0 0
1 −p(x2) −2 +q(x2) Δ2 2Δ Δ2
0 1 −p(x3) −2 +q(x3) 1 +p(x3) 0 Δ2 2Δ Δ2 Δ2 2Δ
0 0 1 −p(x4) −2 +q(x4) 1 +p(x4) Δ2 2Δ Δ2 Δ2 2Δ
000001
α
y0 y1 y2 y3 y4 y5
=
r(x1) r(x2) r(x3) r(x4)
This is now a 6×6 set of linear algebraic system of equations. You can use any of the
y0 y1 y2 y3 y4
y5 11
technique in lectures 5 and 6 to solve for
β
You can slightly reduce the number of unknowns by noting that you do not need to explicitly solve for y0 and y5 because you already know what they are, y0 = α and
y5 = β . So you can rewrite
11p(x1) 20 10p(x1) 0 0 0
Δ2−2Δ −Δ2+q(x1) Δ2+2Δ 0 0
0
0
0
α
0 0 1 −p(x4) −2 +q(x4) 1 +p(x4) Δ2 2Δ Δ2 Δ2 2Δ
y0 r(x )
0
0
0
1 −p(x2) −2 +q(x) 1 +p(x2) 0 Δ2 2Δ Δ2 2 Δ2 2Δ
0 1 −p(x3) −2 +q(x) 1 +p(x3) Δ2 2Δ Δ2 3 Δ2 2Δ
y1 r(x1) y2 =2 y3 r(x3)
y4 r(x4) y5 β
(15.2)
(1 p(x)) yr(x1)−α−1
000001
as
−2 +q(x) 1 +p(x1) 0 Δ2 1Δ22Δ
1 −p(x2) −2 +q(x) 1 +p(x2) Δ2 2Δ Δ2 2 Δ2 2Δ
0 0
1 Δ2 2Δ y2 = r(x2)
0 1 −p(x3) −2 +q(x) 1 +p(x3)
Δ2 0
2Δ Δ2
0
1 −p(x4) −2 +q(x4) Δ2 2Δ Δ2
y3 3 Δ2 2Δ y4
r(x3)
(1 p(x4)) r(x4)−β Δ2 + 2Δ
12
If we define
Then Eq. (15.2) can be written as
Tridiagonal
matrix. Easy to solve see question 3 Assignment 2.
β
The above matrix system was for 5 intervals. In general if you divide a domain up into n intervals,
y1 y2 y3 yn−1 ΔΔΔ Δ
x0 = a x1 x2 x3 xn−1 xn then the system of equations that you need to solve would look like
y0 y1 y2
fi= 1 −p(xi) Δ2 2Δ
gi=− 2 +q(xi) Δ2
hi= 1 +p(xi) Δ2 2Δ
α r(x1)
= r(x2) r (x3) r(x4)
1 0 0 f1 g1 h1
00 00
h20 g3 h3
0 y0 0 y1 0 y2
0 f2 g2 0 0 f3 0 0 0 0 0 0
0 y3 h4 y4
f4 g4
001 y5
y0
10…………0
f1 g1 h1 0 … … 0 0 f2 g2 h2 0 … 0 00⋱⋱⋱⋱0 0 0 ⋱ fn−2 gn−2 hn−2 0 0 0 0 ⋱ fn−1 gn−1 hn−1
0000001
α
r(x1) r (x2)
⋮ r(xn−2)
r (xn−1)
β
13
y⋮ = n−2
yn−1 yn
yn
10…………0 y0 f1 g1 h1 0 … … 0 y1
α r(x1)
r(x2) 00⋱⋱⋱⋱0⋮=⋮
0 f2 g2 h2 0 … 0 y2 0 0 ⋱ fn−2 gn−2 hn−2 0 yn−2
it is (probably) more efficient if you (slightly) reduce the system by putting the value of y0 and yn to the right hand side of the system of equations above.
g1 h1 0 … … y1 r(x1)−f1α f2 g2 h2 0 … y2 r(x2)
0⋱⋱⋱⋱y⋮= ⋮
0 0 0 ⋱ fn−1 gn−1 hn−1 yn−1 0 0 0 0 0 0 1 yn
r(xn−2) r(xn−1) β
0 ⋱ fn−2 gn−2 hn−2 0 0 ⋱ fn−1 gn−1
n−2 r(xn−2) yn−1 r(xn−1) − hn−1β
14
Unequally Spaced Grids.
y0 y1 y2 y3 yn−1
yn = β x0=ax1x2 x3xn−1xn=b
d2y (xi) = yi+1 (xi − xi−1) − yi (xi+1 − xi−1) + yi−1 (xi+1 − xi) d x2 12 (xi+1 − xi−1) (xi+1 − xi) (xi − xi−1)
=α
d2y + p(x) dy + q(x)y = r(x) dx2 dx
y′(xi) = yi+1 − yi−1 xi+1 − xi−1
(15.1)
Sometimes, it is useful to have unequally spaced grids. Using CDA approximation for all the derivatives for unequally spaced data points, Eq (15.1) can be approximated as
with Dirichlet boundary conditions
y0 = α yn = β
(15.3)
− p(xi) + 2 y (x−x )
i+1 i−1 (xi+1 − xi−1) (xi − xi−1) i−1 +(− 2 +q(x))y
(xi+1 − xi) (xi − xi−1) i i
(x−x )
i+1 i−1 (xi+1 − xi−1) (xi+1 − xi) i+1
+ p(xi) + 2 y =r(x)
i
15
If we define
p(x) x −x2x−x ( 2 ) fi=(−xi+1−ixi−1 +(i+1 i−1)(i i−1)) gi= −(xi+1−xi)(xi−xi−1)+q(xi)
h=p(xi)+ 2
(x−x )
i i+1 i−1 (xi+1 − xi−1) (xi+1 − xi) then Eq. (15.3) can be put into matrix form as
10…………0 y0 f1 g1 h1 0 … … 0 y1
α
0 f2 g2 h2 0 … 0 y2 00⋱⋱⋱⋱0⋮=⋮
0 0 ⋱ fn−2 gn−2 hn−2 0 yn−2 0 0 0 ⋱ fn−1 gn−1 hn−1 yn−1
r(xn−2) r(xn−1) β
r(x1)−f1α r(x2)
0 0 0 0 0 0 1 yn and moving the boundary conditions to the RHS gives
r(x1) r(x2)
g1 h1 0 … … y1 f2 g2 h2 0 … y2
0⋱⋱⋱⋱y⋮= ⋮
0 ⋱ fn−2 gn−2 hn−2 n−2 r(xn−2)
0 0 ⋱ fn−1 gn−1 yn−1 r(xn−1) − hn−1β 16
Example L15.1: Solve linear BVP using finite difference. Taken from “Applied Numerical Methods Using MATLAB” by Yang, Cao, Chung and Morris
y′′(x) + 2x y′(x) − x2 y(x) = 0 2
x ∈ [1,2] and the boundary conditions y(1) = 5 and y(2) = 3. Use equally spaced data points. Compare your solution with the exact expression for
y ( x ) = x + x4 2
17
I have just shown that to solve
d2y +p(x)dy +q(x)y=r(x) dx2 dx
y′′(x) + 2x y′(x) − x2 y(x) = 0 2
For this problem, p(x) = 2/x, and q(x) = − 2/x2 The domain is
1≤x≤2
with Dirichlet boundary conditions
y(1)=5 ⇒y0 =5 y(2)=3 ⇒yn =3
and using 2nd order central difference approximation for the first and second derivatives, we need to solve the system below to get the solution for y1,y2,……yn−1.
in the domain
with Dirichlet boundary conditions
a≤x≤b
y(a)=α ⇒y0 =α
y(b)=β ⇒yn =β
and using 2nd order central difference approximation for the first and second derivatives, we need to solve the system below to get the solution for y1,y2,……yn−1.
g1 h1 0 … … y1 f2 g2 h2 0 … y2
r(x1)−f1α r(x2)
⋮
r (xn−2)
r(xn−1) − hn−1β
g1 h1 0 … … y1 f2 g2 h2 0 … y2
r(x1)−5f1 r(x2)
⋮
r (xn−2)
r(xn−1) − 3hn−1
0 ⋱ ⋱ ⋱ ⋱ y⋮ = 0 ⋱ fn−2 gn−2 hn−2 n−2
0 ⋱ ⋱ ⋱ ⋱ y⋮ = 0 ⋱ fn−2 gn−2 hn−2 n−2
0 0 ⋱ fn−1 gn−1 yn−1
0 0 ⋱ fn−1 gn−1 yn−1
18
gi = − Δ2 − x2 2 i2
hi = 1 + 2 Δ2 2Δxi
fi= 1 −p(xi) Δ2 2Δ
g i = − Δ2 + q ( x i ) 2
hi= 1 +p(xi) Δ2 2Δ
fi = 1 − 2 Δ2 2Δxi
Lecture15M01.m
A(1,1)=-(2/Delta.^2)-2./(x(2).^2); A(1,2)=(1/Delta.^2)+(2./x(2))*(1/(2*Delta));
A(i,i-1)=(1/Delta.^2)-(2./x(i+1))*(1/(2*Delta)); A(i,i)=-(2/Delta.^2)-2./(x(i+1).^2); A(i,i+1)=(1/Delta.^2)+(2./x(i+1))*(1/(2*Delta));
A(n-1,n-2)=(1/Delta.^2)-(2./x(n))*(1/(2*Delta)); A(n-1,n-1)=-(2/Delta.^2)-2./(x(n).^2);
C(1)=-alpha*((1/Delta.^2)-(2./x(2))*(1/(2*Delta))); C(n-1)=-beta*((1/Delta.^2)+(2./x(n))*(1/(2*Delta)));
hold off; plot(x,ysol,’bo-‘,’Markersize’,15,’Markerfacecolor’,[0 0 1]); hold on
plot(x,x+4./x.^2,’k-‘,’linewidth’,2);
clear all clear all;
Delta=0.2; x=1:Delta:2; alpha=5.0; beta=3.0;
n=length(x)-1; A=zeros(n-1,n-1); C=zeros(n-1,1); ysol=zeros(size(x));
for i=2:n-2
end
%
% This is a very inefficient way of solving this system of equations
% Matrix [A] has got lots of zeros
% Should use Thomas algorithm and not store all the zeros. %
yinner=A\C;
ysol(1)=alpha; ysol(2:n)=yinner; ysol(n+1)=beta;
xlabel(‘x’,’FontSize’,14) ylabel(‘y(x)’,’FontSize’,14)
22 gi=−Δ2 −x2
i
0 … … y1
1 2 fi=Δ2 −2Δx
1 2 hi=Δ2 +2Δx
g1 h1
i
i
r(x1)−5f1 r(x2)
f2 g2 h2 0 … 0⋱⋱⋱⋱
0 ⋱ fn−2 gn−2 hn−2
0 0
⋱
yn−2 r(xn−2)
n−1 r(xn−1) − 3hn−1
19
fn−1 gn−1
y2 ⋮=⋮
5 4.8 4.6 4.4 4.2 4 3.8 3.6 3.4 3.2 3
Exact solution Computed solution
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 x
y(x)
End of Example L15.1
20
Heated Rod Problem
Solve the boundary value problem
d2T dx2
With boundary conditions
T′(0) = 0 T(1) = 1
0≤x≤1
T=1
=1
Insulated
21
L15.3:
Boundary Value Problem:
Finite Difference Approach
(Mixed (Neumann & Dirichlet) Boundary Conditions)
22
Example L15.2: Solve the following boundary value problem
d 2 T + 1 d T = d ( T 4 − T R4 ) − g ( r ) dr2 r dr
r ∈ [0,R] and the boundary conditions dT/dr(r = 0) = 0 and T(r = R) = TR. Use equally spaced data points.
Neumann boundary condition
Dirichlet boundary condition
23
Heated Disk Problem
A disk is subjected to a radially symmetric heat source, g(r) and is cooled by the conduction to the rim of the disk and radiation to its environment, modelled by the following nonlinear ordinary differential equation
where
and given
d 2 T + 1 d T = d ( T 4 − T R4 ) − g ( r ) dr2 r dr
0≤r≤R T(R) = TR, T′(0) = 0
R=5, d=0.1,
TR =Tb =10,
g(r)=(r−R)2
T
r
24
T0
T1 T2 T3 T4
T5
ΔΔΔΔΔ
r0 = 0 r1 d2T+1dT=d(T4−TR4)−g(r)
d 2 T + 1 d T − d ( T 4 − T R4 ) + ( r − R ) 2 = 0 dr2 r dr
r2 g(r)=(r−R)2
r3 r4 T(R)=T, T′(0)=0
dr2 r dr
R T 5 = T R
r5 = R T1 = T0
Use forward difference approximation. Cannot use central difference approximation because we do not have information at r = r−1
25
dT(r0=0)≈T1−T0 =0 dr Δ
⇒
T′′ r1T1′ T4−T4 (r −R)2 111R10
T ′ ′ 1 T 2′ T 4 − T 4 ( r − R ) 2 2+r2 −d2R+2 =0
T′′ 1T′ T4−T4 (r −R)2
3 r33 3R 3 0
T′′ 1 T4−T4 (r −R)2 4 r4T4′ 4 R 4
ΔΔΔΔΔ
r0 = 0
r1 r2 r3 d2T+1dT−d(T4−TR4)+(r−R)2=0
r4 r5 = R T =T
dr2 r dr
1 0
T′′ r1T1′ T4−T4 (r −R)2 111R10
T ′ ′ 1 T 2′ T 4 − T 4 ( r − R ) 2 2+r2 −d2R+2 =0
T′′ 1T′ T4−T4 (r −R)2
3 r33 3R 3 0
T′′ 1 T4−T4 (r −R)2 4 r4T4′ 4 R 4
T5 = TR
T − 2T + T 1 T′ = 1 T2 − T0
T1′′=2 10r11r12Δd(T4−T4)(r−R)2
Δ2 11T−T1R1
T 2 ′ ′ = T 3 − 2 T 2 + T 1 y 2′ = 3 1 d ( T 24 − T R4 ) ( r 2 − R ) 2
Δ2 r2 r2 2Δ v v
T − 2 T + T v T − T T3′′=4 Δ3 2 1T3′=14 2
(
r − R
( 4 )2
2
r3 r3 2Δ 1 1 y5 − y3
d ( T 34 − T R4 )
d(T4 − T4) 4 R
) 2 r3−R
T ′′ = T5 − 2T4 + T3 4 Δ 2
r 4 T 4′ = r 4 2 Δ
26
T1 T2 T3 T4
T5
T0
ΔΔΔΔΔ
r0 = 0 r1
r2 r3
r4 r5 = R T1 = T0
r1T1′ 44
T′′ 1 11′1R10
T −T (r −R)2 T′′ r T2 T4−T4 (r −R)2
2+2−d2R+2=0 T′′ 1T′ T4−T4 (r −R)2 0
3r333R3 0 T′′ T4−T4 (r −R)2
4 r1T4′ 4 R 4 4
T5 = TR f1(T0,T1,T2,T3,T4,T5)= 2 Δ1 0 +r 22Δ 0 −d(T14−TR4)+(r1−R)2 =0
T 2 − 2 T 1 + T 0 Δ2
1 T 1′ = 1 T 2 − T 0 r1 r1 2Δ
4 4 d(T1 −TR)
( r 1 − R ) 2
T1′′=
T2′′=Δ2 r22r22Δ 2R (3)
1y′= 1 T3−T1 v
d(T4−T4)
(r2−R)2
T3−2T2+T1
T 3 ′ ′ = T 4 − 2 T 3 + T 2 r 1 T 3′ = r 4 2 Δ 2 34 R r − R
Δ2
T 4 ′ ′ = T 5 − 2 T 4 + T 3 Δ2
3
3
2 4 R ( r 4 − R ) 2
4 1 T −T d(T −T )
1 T 4′ = 1 y 5 − y 3
d(T4 −T4)
r4
r4 2Δ
T −2T +T 1 T −T 21
v
v
f2(T0,T1,T2,T3,T4,T5)= T3−2T2+T1 + 1 T3−T1 −d(T24−TR4)+(r2−R)2 =0 Δ2 r22Δ
f3(T0,T1,T2,T3,T4,T5)= T4−2T3+T2 + 1 T4−T2 −d(T34−TR4)+(r3−R)2 =0 Δ2 r32Δ
f4(T0,T1,T2,T3,T4,T5)= T5−2T4+T3 + 1 T5−T3 −d(T4−TR4)+(r4−R)2 =0 Δ2 r42Δ
27
T1 T2 T3 T4
T5
T0
ΔΔΔΔΔ
r0 = 0 r1
r2 r3
r4 r5 = R T1 = T0
This is a set of
4 nonlinear algebraic equations with
4 unknowns, T1, T2, T3, T4. We should be able to solve using Newton-Raphson method or using the fsolve() function in MATLAB.
r1T1′ 44
T′′ 1 11′1R10
T −T (r −R)2 T′′ r T2 T4−T4 (r −R)2
2+2−d2R+2=0 T′′ 1T′ T4−T4 (r −R)2 0
3r333R3 0 T′′ T4−T4 (r −R)2
4 r1T4′ 4 R 4 4
T5 = TR f1(T1,T2,T3,T4)= 2 Δ1 1 +r 22Δ 1 −d(T14−TR4)+(r1−R)2 =0
T 2 − 2 T 1 + T 0 Δ2
1 T 1′ = 1 T 2 − T 0 r1 r1 2Δ
4 4 d(T1 −TR)
( r 1 − R ) 2
T1′′=
T2′′=Δ2 r22r22Δ 2R (3)
1y′= 1 T3−T1 v
d(T4−T4)
(r2−R)2
T3−2T2+T1
T 3 ′ ′ = T 4 − 2 T 3 + T 2 r 1 T 3′ = r 4 2 Δ 2 34 R r − R
Δ2
T 4 ′ ′ = T 5 − 2 T 4 + T 3 Δ2
3
3
2 R ( r 4 − R ) 2
r4
r4 2Δ
T −2T +T 1 T −T 21
4 1 T −T d(T −T )
1 T 4′ = 1 y 5 − y 3
d(T4 −T4) 4
v
v
f2(T1,T2,T3,T4)= T3−2T2+T1 + 1 T3−T1 −d(T24−TR4)+(r2−R)2 =0 Δ2 r22Δ
f3(T1,T2,T3,T4)= T4−2T3+T2 + 1 T4−T2 −d(T34−TR4)+(r3−R)2 =0 Δ2 r32Δ
f4(T1,T2,T3,T4)= TR−2T4+T3 + 1 TR−T3 −d(T4−TR4)+(r4−R)2 =0 Δ2 r42Δ
28
T1 T2 T3 T4
T5
T0
If we now have n + 1 (and n intervals)
ΔΔΔ Δ
Tn
Tn−1
T0
T2 T3
T1
r0 = 0 r1
r2 r3
d 2 T + 1 d T = d ( T 4 − T R4 ) − g ( r )
rn−1 rn = R
: r = r1
: r = r2,r3,…rn−2
: r = rn−1
T′(0) = 0 T1 = T0
dr2 r dr
T(R) = TR Tn = TR
f 1 = T 2 − 2 T 1 + T 1 + 1 T 2 − T 1 − d ( T 14 − T R4 ) + ( r 1 − R ) 2 = 0 Δ2 r12Δ
f i = T i + 1 − 2 T i + T i − 1 + 1 T i + 1 − T i − 1 − d ( T i 4 − T R4 ) + ( r i − R ) 2 = 0 Δ2 ri2Δ
fn−1=TR−2Tn−1+Tn−2 + 1 TR−Tn−2 −d(T4 −TR4)+(rn−1−R)2=0 Δ2 rn−1 2Δ n−1
29
Lecture15M02.m
clear all d=0.1;TR=10;R=5;
R=5, d=0.1, n=length(r)-1; %n is the number of intervals
T=zeros(1,n+1);
%Note that I am only solving for the inner variable!
TR =10
Delta=0.2; r=0:Delta:R;
Remember that we only need
to solve for the inner variables since from the boundary conditions we know thatT0 =T1 andTn =TR =10
Tn
rn = R
Tinner = fsolve(@(x) HeatedDiskExample(x,r(2:n),R,d,TR,Delta),10*ones(1,n-1));
%Putting the boundary values back in the solution vector.
T(1)=Tinner(1); T(2:n)=Tinner(1:end); T(n+1)=TR;
plot(r,T,’ko-‘,’linewidth’,2,’Markersize’,10,’Markerfacecolor’,[0 0 1]) xlabel(‘r’,’FontSize’,14)
ΔΔΔ Δ
r2 r3 rn−1
r0 = 0 r1 f(1) = (x(2)-2*x(1)+x(1))/Delta^2+(1./(2*Delta*r(1)))*(x(2)-x(1))-d*(x(1).^4-TR.^4)+(r(1)-R).^2;
ylabel(‘T’,’FontSize’,14)
function f = HeatedDiskExample(x,r,R,d,TR,Delta) for i=2:length(r)-1
end
k=length(r);
f(k) = (TR-2*x(k)+x(k-1))/Delta^2+(1./(2*Delta*r(k)))*(TR-x(k-1))-d*(x(k).^4-TR.^4)+(r(k)-R).^2;
f(i) = (x(i+1)-2*x(i)+x(i-1))/Delta^2+(1./(2*Delta*r(i)))*(x(i+1)-x(i-1))-d*(x(i).^4-TR.^4)+(r(i)-R).^2;
end
f1 = T2−2T1+T1 + 1 T2−T1 −d(T14−TR4)+(r1−R)2 =0 Δ2 r12Δ
n−1 TR − 2Tn−1 + Tn−2 1 TR − Tn−2 n−1 R ( n−1 )2
f = Δ2 +rn−1 2Δ −d(T4 −T4)+r −R =0
T0
T1
f i = T i + 1 − 2 T i + T i − 1 + 1 T i + 1 − T i − 1 − d ( T i 4 − T R4 ) + ( r i − R ) 2 = 0 Δ2 ri2Δ
30
T2 T3
Tn−1
Output
10.06
10.05
10.04
10.03
10.02
10.01
10
31
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
r
T
End of Example L15.2
32