392 Boundary value problems
37 k2 = k2 + dk;
38 else
39 fprintf(‘Amazingly got r = 1 to numerical …
precision.\n’)
50 n = 101;
51 dx = 1/(n-1);
52 x = linspace(0,1,n);
53 A = zeros(2*n);
54 b = zeros(2*n,1);
65 A(2*i-1,2*i-3) = 1;
66 A(2*i-1,2*i-1) = -2 – dx*dx*k1/Da;
67 A(2*i-1,2*i) = dx*dx*k2/Da;
68 A(2*i-1,2*i+1) = 1;
69 %equation for B
80
81 %right boundary condition for B
82 A(2*n,2*n-2) = 2;
83 A(2*n,2*n-1) = dx*dx*k1/Db;
84 A(2*n,2*n) = -2 – dx*dx*k2/Db;
Boundary value problems 393
The text output is:
1k2 = 0.000000 r = -0.000000
2
3****************************************
14 k2 = 9.000000 r = 0.8572269
15 k2 = 10.000000 r = 0.9553184
16 k2 = 11.000000 r = 1.0537708
17 Went too far. Reduce step size to dk = 1.0e-01.
18
29 k2 = 10.420000 r = 0.9966269
30 k2 = 10.430000 r = 0.9976112
31 k2 = 10.440000 r = 0.9985955
32 k2 = 10.450000 r = 0.9995798
33 k2 = 10.460000 r = 1.0005642
44 k2 = 10.454000 r = 0.9999736
45 k2 = 10.454100 r = 0.9999834
46 k2 = 10.454200 r = 0.9999933
47 k2 = 10.454300 r = 1.0000031
48 Went too far. Reduce step size to dk = 1.0e-05.
58 Converged to within tolerance of 1.0e-06.
59
60 ****************************************
(c) The Matlab script is:
1function s12c11p123c
12 Dbplot = zeros(npts,1);
13 k2plot = zeros(npts,1);
14
15 %get the value of k2 where the wall concentrations are …
equal for a given
26 ylabel(‘$k 2ˆ*$’,‘FontSize’,14)
27 saveas(h,‘s12c11p123 solution figure2.eps’,‘psc2’)
28
29
30
Boundary value problems 395
43 fprintf(‘k2 = %8.6f \t r = %8.7f \n’,k2,r)
%2.1e.\n\n’,dk)
50 elseif r<1
51 k2 = k2 + dk;
52 else
53 fprintf(‘Amazingly got r = 1 to numerical …
63 n = 101;
64 dx = 1/(n-1);
65 x = linspace(0,1,n);
66 A = zeros(2*n);
67 b = zeros(2*n,1);
78 A(2*i-1,2*i-3) = 1;
79 A(2*i-1,2*i-1) = -2 – dx*dx*k1/Da;
80 A(2*i-1,2*i) = dx*dx*k2/Da;
81 A(2*i-1,2*i+1) = 1;
82 %equation for B
93
94 %right boundary condition for B
95 A(2*n,2*n-2) = 2;
96 A(2*n,2*n-1) = dx*dx*k1/Db;
396 Boundary value problems
Boundary value problems 397
(b) For a given size L, we need n=L/∆z+1 nodes. For node i=1, the boundary
condition gives w1=75. For the interior nodes i=2 : n−1, we have
(∆z)2!+ziwi+1−wi−1
2∆z−z2
The right boundary condition gives wn+1=wn−1, so the equation we need to
(∆z)2!−z2
The Matlab script is:
1function s13c10p123b
2clc
3close all
4set(0,‘defaulttextinterpreter’,‘latex’)
5
6
7L = 0.25; %25 cm long
8t = 0.05; %5 cm thickness
19 A = zeros(npts);
20 b = zeros(npts,1);
21
22 %write the first row
23 A(1,1) = 1; b(1) = 75; %initial temperature 100 – 25 (ref)
34 %write the last row
35 z=z+dz;%location of the node
36 A(npts,npts-1) = 2*z/dzˆ2;
37 A(npts,npts) = -2*z/dzˆ2 – zˆ2;
38
48 ylabel(‘Temperature (C)’,‘FontSize’,14)
49 title(‘Solution to s13c10p123b’,‘FontSize’,14)
50 saveas(h,‘s13c10p123 solution figure1.eps’,‘psc2’)
The output is:
0 5 10 15 20 25
65
70
75
80
85
90
95
100
Dis t anc e (c m )
Te mp e r at ur e ( C )
Solut ion t o s 13c 10p123b
(c) The Matlab script is:
1function s13c10p123c
2clc
3close all
14 %quick start to find an estimate for T
15 fprintf(‘Making a scan with dL = %6.4f m\n’,dL)
16 fprintf(‘================================\n’)
17 fprintf(‘L = %6.4f \t Tend = %6.5f \n’,L,Tend)
18 while Tend >Tinfinity +1
Boundary value problems 399
26 fprintf(‘%6.4f m\n\n’,L)
27
28
29 %decrease dL by a lot
40
41
42
43 function out = getT(L)
44 t = 0.05; %5 cm thickness
54 %initialize the variables
55 A = zeros(npts);
56 b = zeros(npts,1);
57
58 %write the first row
69
70 %write the last row
71 z=z+dz;%location of the node
72 A(npts,npts-1) = 2*z/dzˆ2;
73 A(npts,npts) = -2*z/dzˆ2 – zˆ2;
400 Boundary value problems
1Making a scan with dL = 0.0500 m
2================================
3L = 0.0500 Tend = 98.25493
14 L = 0.6000 Tend = 35.56995
15 L = 0.6500 Tend = 33.90671
16 L = 0.7000 Tend = 32.55023
17 L = 0.7500 Tend = 31.43405
18 L = 0.8000 Tend = 30.50903
29 L = 1.3500 Tend = 26.22935
30 L = 1.4000 Tend = 26.08571
31 L = 1.4500 Tend = 25.96020
32
33
44 L = 1.4060 Tend = 26.06984
45 L = 1.4070 Tend = 26.06719
46 L = 1.4080 Tend = 26.06453
47 L = 1.4090 Tend = 26.06189
48 L = 1.4100 Tend = 26.05925
Boundary value problems 401
50 L = 1.4120 Tend = 26.05411
51 L = 1.4130 Tend = 26.05150
52 L = 1.4140 Tend = 26.04889
53 L = 1.4150 Tend = 26.04629
7Partial differential equations
414 Partial differential equations
Problems
(7.1)
(7.2) The total number of nodes is n2. The interior nodes are (n−2)(n−2) in size and have
(7.4) The IVP for the left node i=1 requires that we solve for the fictitious node from
dc
dx=0
which becomes
c2−c0
The problem is linear, so you can use any of the linear solvers (Gauss elimination,
LU decomposition, relaxation methods) to compute the value at the next time step.
(7.5) Discretize the interior nodes with centered finite differences
dt=0; c(0)
We need to use a fictitious node at the right boundary. With the derivative equal to
zero, we get cn+1−cn−1/(2∆x)=0, which means cn+1=cn−1. The equation at the
right boundary is then
dt=2D
(∆x)2(cn−1−cn)−kc2
Partial differential equations 415
If we integrate with forward Euler, then the left node is
The interior nodes are
(∆x)2(c(k)
and the right node is
(∆x)2(c(k)
If we integrate with backward Euler, the residual at the first node is
(∆x)2c(k+1)
(∆x)2+1!+hk(c(k+1)
The Jacobian elements for the first row are J11 =1 and zero otherwise. For interior
nodes, there are three non-zero entries:
Ji,i−1=−hD
(∆x)2+2hkc(k+1)
(7.6) For the interior nodes, the equation is
dci
dt=D
(∆x)2(ci+1−2ci+ci−1)
The normal vector points outwards, so n=−ixon the left boundary and n=ixon the
right boundary. So the left boundary condition is
(7.7) From the fictitious node equation on the left, we have
c2−c0
2∆x=0
which gives c0=c2. The system of equations is thus
dc1
dt=c3 1
∆x2!−c2 2
∆x2+k!+c1 1
∆x2!
dc3
dt=0
Partial differential equations 417
If we further substitute the values for the parameters, ∆x=1/2 so 1/∆x2=4,
4(0)
Eliminating c3from the second equation with c3=1, we have the system
7
2c1−2c2=1
−c1+7
(7.8) The ODEs you need to solve in the interior (2-5) are
dci
dt =Dci+1−2ci+ci−1
(∆x)2−kci
418 Partial differential equations
The left node has the equation
t=0.1 1.4 1.4 1.4 1.4 -8.6 1
(7.9) If we start numbering the 9 nodes at the lower-left corner, then the interior node
equation for c5is
c6−2c5+c4
(∆x)2+c8−2c5+c2
(∆y)2=0
If we use a square grid with ∆x= ∆y, then
Since the system is symmetric, we still get c5=3/2. In contrast, if the upper boundary
was c=3, then we would break the symmetry and get
If we go back to the square grid ǫ=1 and replace the c=3 condition at the top
with no flux, the Dirchelet condition prevails on the corneres so c7=1 and c9=2.
For the other corners, we take the average of the boundary temperatures, c1=1 and
c3=3/2. We need to define an auxiliary node c11. The boundary condition on node
8 is n· ∇c=0. The normal vector points in yso we have ∂c/∂y=0. This gives the
result c11=c5. Using the differential equation for node 8 and ∆x= ∆ywe get
(7.10) 5/(nxny), nx+ny+1 and nxny
(7.11) The residual for node 7 is
c8−4c7+c6+c11 +c3=0
and the non-zero entries in the Jacobian are
J7,8=J7,6=J7,11 =J7,3=1
420 Partial differential equations
and the non-zero entry in the Jacobian is
J14,14 =1
(7.12) (a) Finite differences with a banded solver for the linear problem. The matrix size
is (2n)2with a bandwidth of n. So the time is n4. If we increase the nodes by a
factor of 10, the time goes up by 104and takes 103seconds.
(b) Method of lines (with finite differences) and an explicit solver like explicit Euler.
For explicit Euler, the maximum step size is
If we use a step size of 10−6, we need 107steps to reach t=10. If each step takes
10−2seconds, then the total time is at least 105seconds.
(7.13) The method of lines gives
The 4th variable z4=T2. So we need to write out the equation
From implicit Euler, we would have
There are 4 non-zero entries in the Jacobian:
∆x2+φγe−β/z4,i+1
z2
4,i+1
J4,6=−h
∆x2
(7.14) To get a higher-order method, you will need to use 5 pts centered around the node
and expand up to the fourth derivative,
Partial differential equations 421
For a centered finite differences, we need four expansions:
We want to get cancellation of all of the even ordered derivatives, which we can do
by forming the sum
S=α[f(xi+1)−f(xi−1)] +β[f(xi+2)−f(xi−2)]
WIth some obvious simplifications, we get
Computer Problems
(7.15) The files for this problem are in the folder s13c11p3 matlab.
The method of lines equations with finite differences are
with c1=1 and cn=0. Since these values do not change, their derivatives need to be
set to zero in the method of lines.
422 Partial differential equations
1function s13c11p3
2clc
3close all
12 dt = 0.00001;
13 tmax = 0.1;
23 for j = tstart:tend
24 c = RK4(c,dt,dx,npts); %integrate
25 end
26 cplot(i,:) = c; %store output
27 end
38 function out = RK4(y,h,dx,npts)
39 k1 = feval(y,dx,npts);
40 k2 = feval(y+h*k1/2,dx,npts);
41 k3 = feval(y+h*k2/2,dx,npts);
42 k4 = feval(y+h*k3,dx,npts);