[Content_Types].xml
_rels/.rels
matlab/document.xml
matlab/output.xml
metadata/coreProperties.xml
metadata/mwcoreProperties.xml
metadata/mwcorePropertiesExtension.xml
metadata/mwcorePropertiesReleaseInfo.xml
The Finite Difference Method for Solving Boundary Value Problems In this livescript, you will learn how To solve boundary value problems using finite difference methods. We’ll continue with our example on steady one dimensional heat diffusion \frac{d^{2}T}{dx^{2}}=0 with the boundary conditions T(0)=0 and T(1)=1 . The process of solving any boundary value problem using the finite difference method boils down to a couple of simple steps Generate a mesh with a sufficient number of nodes. Discretise the ODE using a finite difference scheme. Construct the system of equations. Implement the boundary conditions. Solve for the unknown nodes. Step 1. Mesh generation. To generate an equally spaced mesh, we can just use \texttt{linspace} . N = 10;
xmesh = linspace(0,1,N+1) However be aware that in many physical problems, the physics of the problem will determine how we should distribute our points to most efficiently solve the problem. Step 2. Discretise the ODE. It’s generally a safe bet to use a central difference scheme to discretise diffusion type problems. However, be wary when convective effects are a part of the problem, as central differencing can cause unphysical oscillations to appear in the solution. So applying central differencing to our ODE gives \frac{T_{i-1}-2T_{i}+T_{i+1}}{\Delta^{2}}=0 for i=2,\ldots,N . NOTE: In the lecture notes the index starts from 0 . However, due to the indexing in MATLAB, we’ll shift it up by 1 . We can equivalently write this as T_{i-1}-2T_{i}+T_{i+1}=0 In matrix form, this looks like \pmatrix{
\ddots& & & \vdots && & \cdot^{\cdot^{\cdot}}\cr
\cdots& 0 & 1 & -2 & 1 & 0 & \cdots \cr
\cdot^{\cdot^{\cdot}} & & & \vdots & & & \ddots
}
\pmatrix{
\vdots \cr T_{i-1} \cr T_{i} \cr T_{i+1} \cr \vdots
}
=
\pmatrix{
\vdots \cr
0 \cr
0 \cr
0 \cr
\vdots
} Step 3. System Assembly. The first step as always is to initialise any arrays. C = zeros(N-1,N+1);
b = zeros(N-1,1); Then, the interior nodes of the mesh can be constructed by using the following loop. for k = 1:N-1
C(k,k:k+2) = [1,-2,1];
b(k) = 0;
end Step 4. Boundary Conditions. Since the boundary conditions are by definition set, we can write them as T_{0}=0\\T_{N}=1 To apply these, we’ll expand the matrix to give T_{0}
\pmatrix{
1 \cr
0 \cr
\vdots \cr
0 \cr
0 \cr
}
+
\pmatrix{
-2 & 1 & & &\cr
1 & -2 & 1 & & \cr
& \ddots & \ddots & \ddots & \cr
& & 1 & -2 & 1 \cr
& & & 1 & -2
}
\pmatrix{
T_{1}\cr
T_{2}\cr
\vdots\cr
T_{N-2}\cr
T_{N-1}\cr
}
+
T_{N}
\pmatrix{
0 \cr
0 \cr
\vdots \cr
0 \cr
1 \cr
}
=
\pmatrix{
b_{1} \cr
b_{2} \cr
\vdots \cr
b_{N-2} \cr
b_{N-1} \cr
} Rearranging for the unknowns gives the linear system \pmatrix{
-2 & 1 & & &\cr
1 & -2 & 1 & & \cr
& \ddots & \ddots & \ddots & \cr
& & 1 & -2 & 1 \cr
& & & 1 & -2
}
\pmatrix{
T_{1}\cr
T_{2}\cr
\vdots\cr
T_{N-2}\cr
T_{N-1}\cr
}
=
\pmatrix{
b_{1}-T_{0} \cr
b_{2} \cr
\vdots \cr
b_{N-2} \cr
b_{N-1}-T_{N} \cr
} And a piece of code that adds this to our [C] matrix and \{b\} vector is C_sys = C(:,2:end-1);
b(1) = b(1)-0;
b(end) = b(end)-1; Step 5. Solving the System. Before we just use \texttt{\\} to solve our system, let’s just have a look at its structure to see what we’re dealing with C_sys We can see that [C] has a tridiagonal structure, and thus may be efficiently solved using the Thomas algorithm . If we were to increase the order of the central difference method to, say, 4th order \frac{-T_{i-2}+16T_{1}-30T_{i}+16T_{i+1}-T_{i+2}}{12\Delta^{2}}=0 our [C] matrix would just be pentadiagonal, which can still be solved efficiently by modifying the Thomas algorithm. For simplicity sake we’ll just continue with \texttt{\\} , which gives the solution. T_sys = C_sys\b
T = [0;T_sys;1];
plot(xmesh,T)
xlabel(‘x’)
ylabel(‘Temperature’) (a) Modify the code to use the Thomas algorithm. See which piece of code performs better.
manual code error 0.4 true false 0 14 14 false true 1 15 15 true false 2 31 31 false false 3 32 32 false true 4 34 37 true false 5 47 47 false false 6 48 48 false true 7 49 49 true false 8 53 53 false false 9 59 59 false false 10 60 60 false false 11 61 61 false false 12 62 62 false true 13 63 63
2020-04-30T15:12:13Z 2020-09-27T02:53:23Z
application/vnd.mathworks.matlab.code MATLAB Code R2020b
9.9.0.1444674 64b92d2a-6b34-4917-aca3-62a6fc59ebd9
9.9.0.1467703
R2020b
Aug 26 2020
2314982500