### Related Posts

The ODE system above cannot be used for $$u_{0}^{\prime}$$ since that equation involves some quantity $$u_{-1}^{\prime}$$ outside the domain. We remark that the temperature in a fluid is influenced not only by diffusion, but also by the flow of the liquid. The surface temperature at the ground shows daily and seasonal oscillations. 0 & 0 & 0 & 0 & \dots & 0 & 1 & -2 & 1 & 0 \\ pp 153-175 | where $$\alpha$$ is the thermal conductivity of the rod and $$\sigma (x,t)$$ is a heat source present along the rod. # notice how we multiply numpy arrays pointwise. The surface corresponds to x = 0 and the x axis point downwards into the ground. A solution to a boundary value problem is a solution to the differential equation which also satisfies the boundary conditions. The heat equation around grid node $$1$$ is then modified as: The effect of the Neumann boundary condition is two-fold: it modifies the left-hand side matrix coefficients and the right-hand side source term. The type and number of such conditions depend on the type of equation. To avoid oscillations one must have $$\Delta t$$ at maximum twice the stability limit of the Forward Euler method. In an introductory book like this, nowhere near full justice to the subject can be made. Actually, this reduces the work from the order N 3 to the order N. In one-dimensional diffusion problems, the savings of using a tridiagonal matrix are modest in practice, since the matrices are very small anyway. What about the source term g in our example with temperature distribution in a rod? The dsolve function finds a value of C1 that satisfies the condition. Trying out some simple ones first, like, The simplest implicit method is the Backward Euler scheme, which puts no restrictions on, $$\frac{u^{n+1}-u^{n}}{\Delta t}=f(u^{n+1},t_{n+1})\thinspace.$$, In our case, we have a system of linear ODEs (, $$\frac{u_{0}^{n+1}-u_{0}^{n}}{\Delta t} =s^{\prime}(t_{n+1}),$$, $$\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t} =\frac{\beta}{\Delta x^{2}}(u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1})+g_{i}(t_{n+1}),$$, $$\frac{u_{N}^{n+1}-u_{N}^{n}}{\Delta t} =\frac{2\beta}{\Delta x^{2}}(u_{N-1}^{n+1}-u_{N}^{n+1})+g_{i}(t_{n+1})\thinspace.$$, $$u_{0}^{n+1} =u_{0}^{n}+\Delta t\,s^{\prime}(t_{n+1}),$$, $$u_{1}^{n+1}-\Delta t\frac{\beta}{\Delta x^{2}}(u_{2}^{n+1}-2u_{1}^{n+1}+u_{0}^{n+1}) =u_{1}^{n}+\Delta t\,g_{1}(t_{n+1}),$$, $$u_{2}^{n+1}-\Delta t\frac{2\beta}{\Delta x^{2}}(u_{1}^{n+1}-u_{2}^{n+1}) =u_{2}^{n}+\Delta t\,g_{2}(t_{n+1})\thinspace.$$, A system of linear equations like this, is usually written on matrix form, $$A=\left(\begin{array}[]{ccc}1&0&0\\ -\Delta t\frac{\beta}{\Delta x^{2}}&1+2\Delta t\frac{\beta}{\Delta x^{2}}&-\Delta t\frac{\beta}{\Delta x^{2}}\\ 0&-\Delta t\frac{2\beta}{\Delta x^{2}}&1+\Delta t\frac{2\beta}{\Delta x^{2}}\end{array}\right)$$, $$A_{i,i-1} =-\Delta t\frac{\beta}{\Delta x^{2}},\quad i=2,\ldots,N-1$$, $$A_{i,i+1} =-\Delta t\frac{\beta}{\Delta x^{2}},\quad i=2,\ldots,N-1$$, $$A_{i,i} =1+2\Delta t\frac{\beta}{\Delta x^{2}},\quad i=2,\ldots,N-1$$, $$A_{N,N-1} =-\Delta t\frac{2\beta}{\Delta x^{2}}$$, $$A_{N,N} =1+\Delta t\frac{2\beta}{\Delta x^{2}}$$, If we want to apply general methods for systems of ODEs on the form, $$K_{i,i-1} =\frac{\beta}{\Delta x^{2}},\quad i=2,\ldots,N-1$$, $$K_{i,i+1} =\frac{\beta}{\Delta x^{2}},\quad i=2,\ldots,N-1$$, $$K_{i,i} =-\frac{2\beta}{\Delta x^{2}},\quad i=2,\ldots,N-1$$, $$K_{N,N-1} =\frac{2\beta}{\Delta x^{2}}$$, $$K_{N,N} =-\frac{2\beta}{\Delta x^{2}}$$, $$u(0,t)=T_{0}+T_{a}\sin\left(\frac{2\pi}{P}t\right),$$, Show that the present problem has an analytical solution of the form, An equally stable, but more accurate method than the Backward Euler scheme, is the so-called 2-step backward scheme, which for an ODE, $$\frac{3u^{n+1}-4u^{n}+u^{n-1}}{2\Delta t}=f(u^{n+1},t_{n+1})\thinspace.$$, We consider the same problem as in Exercise, $$E=\sqrt{\Delta x\Delta t\sum_{i}\sum_{n}(U_{i}^{n}-u_{i}^{n})^{2}}\thinspace.$$, The Crank-Nicolson method for ODEs is very popular when combined with diffusion equations. 'Heat equation - Homogeneous Dirichlet boundary conditions'. In such a case, we can split the domain in two and compute u in only one half, $$[-1,0]$$ or $$[0,1]$$. Dirichlet boundary conditions result in the modification of the right-hand side of the equation, while Neumann boundary conditions result into the modification of both the left-hand side and the right-side of the equation. We therefore have a boundary condition $$u(0,t)=s(t)$$. Rather, one must resort to more efficient storage formats and algorithms tailored to such formats, but this is beyond the scope of the present text. The time interval for simulation and the time step depend crucially on the values for β and L, which can vary significantly from case to case. Commonly used boundary conditions are. b_{nx-3} \\ Solving Partial Differential Equations In a partial differential equation (PDE), the function being solved for depends on several variables, and the differential equation can include partial derivatives taken with respect to each of the variables. # is needed. Apply the Crank-Nicolson method in time to the ODE system for a one-dimensional diffusion equation. 3. 0 & 0 & 0 & 0 & \dots & 1 & -2 & 1 & 0 & 0 \\ u (x, 0) = sin (π x). Mathematically, (with the temperature in Kelvin) this example has $$I(x)=283$$ K, except at the end point: $$I(0)=323$$ K, $$s(t)=323$$ K, and g = 0. For convenience, we start with importing some modules needed below: Let’s consider a rod made of heat conducting material. \begin{pmatrix} T_{nx-3} \\ This is nothing but a system of ordinary differential equations in $$N-1$$ unknowns $$u_{1}(t),\ldots,u_{N-1}(t)$$! This service is more advanced with JavaScript available, Programming for Computations - MATLAB/Octave However, there are occasions when you need to take larger time steps with the diffusion equation, especially if interest is in the long-term behavior as $$t\rightarrow\infty$$. However, partial differential equations constitute a non-trivial topic where mathematical and programming mistakes come easy. In two- and three-dimensional PDE problems, however, one cannot afford dense square matrices. Physically this corresponds to specifying the heat flux entering or exiting the rod at the boundaries. One could think of chemical reactions at a microscopic level in some materials as a reason to include g. However, in most applications with temperature evolution, g is zero and heat generation usually takes place at the boundary (as in our example with $$u(0,t)=s(t)$$). What happens inside the rod? The subject of partial differential equations (PDEs) is enormous. Using a Forward Euler scheme with small time steps is typically inappropriate in such situations because the solution changes more and more slowly, but the time step must still be kept small, and it takes ‘‘forever’’ to approach the stationary state. Looking at the entries of the K matrix, we realize that there are at maximum three entries different from zero in each row. \begin{pmatrix} u (0, t) = 0, π e-t + ∂ u ∂ x (1, t) = 0. We can run it with any $$\Delta t$$ we want, its size just impacts the accuracy of the first steps. Reformulate the problem in Exercise 5.6 such that we compute only for $$x\in[0,1]$$. The system can be solved by inverting $$\tilde A_{ij}$$ to get: Inverting matrices numerically is time consuming for large-size matrices. Finally, u(i) has the same indices as rhs: u(2:N). In the previous solution, the constant C1 appears because no condition was specified. 4.2.6. The relevant object name is ThetaRule: Consider the physical application from Sect. Let (5.38) be valid at mesh points x i in space, discretize $$u^{\prime\prime}$$ by a finite difference, and set up a system of equations for the point values u i ,$$i=0,\ldots,N$$, where u i is the approximation at mesh point x i . We have to introduce a discrete version of the condition $$T'(0)=2$$. The initial condition $$u(x,0)=I(x)$$ translates to an initial condition for every unknown function $$u_{i}(t)$$: $$u_{i}(0)=I(x_{i})$$, $$i=0,\ldots,N$$. \begin{pmatrix} The boundary condition reads $$u(0,t)=s(t)$$. A nice feature with having a problem defined as a system of ODEs is that we have a rich set of numerical methods available. Zhi‐Zhong Sun, Weizhong Dai, A new higher‐order accurate numerical method for solving heat conduction in a double‐layered film with the neumann boundary condition, Numerical Methods for Partial Differential Equations, 10.1002/num.21870, 30, 4, (1291-1314), (2014). It would be much more efficient to store the matrix as a tridiagonal matrix and apply a specialized Gaussian elimination solver for tridiagonal systems. This operation is performed with the help of the numpy.dot function that allows many sorts of vector and matrix multiplications. In the previous notebook we have defined $$A_{ij}$$ for the centered second-order accurate second-order derivative as: Let’s see how to modify this matrix to take into account the boundary conditions. One dimensional heat equation: implicit methods Iterative methods 1. Indeed, in many problems, the loss of accuracy used for the boundary conditions would degrade the accuracy of the solution throughout the domain. Such situations can be dealt with if we have measurements of u, but the mathematical framework is much more complicated. In a later chapter of this course we will explain how to obtain approximate inverses for large systems. Also note that the rhs function relies on access to global variables beta, dx, L, and x, and global functions dsdt, g, and dudx. One such class is partial differential equations (PDEs). $$\varrho=2.7\cdot 10^{3}\hbox{ kg/m}^{3}$$, $$\kappa=200\,\,\frac{\hbox{W}}{\hbox{mK}}$$, $$\beta=\kappa/(\varrho c)=8.2\cdot 10^{-5}\hbox{ m}^{2}/\hbox{s}$$, Exercise 5.1 (Simulate a diffusion equation by hand), Exercise 5.2 (Compute temperature variations in the ground), Exercise 5.4 (Explore adaptive and implicit methods), Exercise 5.6 (Compute the diffusion of a Gaussian peak), Exercise 5.7 (Vectorize a function for computing the area of a polygon), $$x_{1}y_{2}+x_{2}y_{3}+\cdots+x_{n-1}y_{n}=\sum_{i=0}^{n-1}x_{i}y_{i+1}$$, Exercise 5.10 (Solve a two-point boundary value problem), http://creativecommons.org/licenses/by-nc/4.0/, Department of Process, Energy and Environmental Technology, https://doi.org/10.1007/978-3-319-32452-4_5, Texts in Computational Science and Engineering. We also have briefly discussed the usage of two functions from scipy and numpy to respectively invert matrices and perform array multiplications. Odespy requires the problem to be formulated in Python code. Vectorize the implementation of the function for computing the area of a polygon in Exercise  2.6. Without them, the solution is not unique, and no numerical method will work. 0 & 0 & 0 & 0 & \dots & 0 & -1 & 4 & -5 & 2 You may use the Forward Euler method in time. Notice that the formula $$x_{1}y_{2}+x_{2}y_{3}+\cdots+x_{n-1}y_{n}=\sum_{i=0}^{n-1}x_{i}y_{i+1}$$ is the dot product of two vectors, x(1:end-1) and y(2,end), which can be computed as dot(x(1:end-1), y(2,end)), or more explicitly as sum(x(1:end-1). Similarly, u(i-1) corresponds to all inner u values displaced one index to the left: u(1:N-1). Partial Differential Equations Version 11 adds extensive support for symbolic solutions of boundary value problems related to classical and modern PDEs. If we look back at equation (31), we have in full generality: If we set $$T_0=1$$, this equation becomes: We observe that compared to our previous setup, the left-hand side has not changed. system of ordinary differential equations. 4.3.6). You should have a look at its documentation page. One dimensional heat equation 4. 0 & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & 0 \\ As the loop index i runs from 2 to N, the u(i+1) term will cover all the inner u values displaced one index to the right (compared to 2:N), i.e., u(3:N+1). Solving Differential Equations online. For partial di erential equations (PDEs), we need to know the initial values and extra information about the behaviour of the solution u(x;t) at the boundary of the spatial domain (i.e. For this particular equation we also need to make sure the initial condition is $$u_{0}(0)=s(0)$$ (otherwise nothing will happen: we get u = 283 K forever). We can compare it with the exact solution $$T(x)=\displaystyle\frac{1}{2}x(1-x)$$, which obviously satisfies the required boundary conditions. Run this case with the θ rule and $$\theta=1/2$$ for the following values of $$\Delta t$$: 0.001, 0.01, 0.05. We should also mention that the diffusion equation may appear after simplifying more complicated partial differential equations. The partial differential equations to be discussed include •parabolic equations, •elliptic equations, •hyperbolic conservation laws. To close this equation, some boundary conditions at both ends of the rod need to be specified. Nonlinear partial differential equation with random Neumann boundary conditions are considered. # The first line of A needs to be modified for the Neumann boundary condition. 0 & 0 & 1 & -2 & 1 & \dots & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots & 0 & 0 & 0 & 1 & -2 The Backward Euler method with $$\Delta t=0.001$$, The backward 2-step method with $$\Delta t=0.001$$, The backward 2-step method with $$\Delta t=0.01$$. For such applications, the equation is known as the heat equation. 0 & 0 & 1 & -2 & 1 & \dots & 0 & 0 & 0 & 0 \\ In this first example, we apply homogeneous Dirichlet boundary conditions at both ends of the domain (i.e. Plot both the numerical and analytical solution. It turns out that solutions, $$3(x-L)=0+g(x,t)\quad\Rightarrow\quad g(x,t)=3(x-L)\thinspace.$$, First we need to generalize our method to handle, $$\frac{u_{N+1}(t)-u_{N-1}(t)}{2\Delta x}=\gamma\quad\Rightarrow\quad u_{N+1}=u_{N-1}+2\gamma\Delta x,$$, $$\frac{du_{N}(t)}{dt}=\beta\frac{2u_{N-1}(t)+2\gamma\Delta x-2u_{N}(t)}{\Delta x^{2}}+g_{N}(t)\thinspace.$$. This is to be expected, as we now have $$nx-2$$ unknowns. # Computation of the solution using numpy.dot (boundary nodes not included). At the surface, the temperature has then fallen. Implement the θ rule with aid of the Odespy package. b_1 \\ A better start is therefore to address a carefully designed test example where we can check that the method works. The type and number of such conditions depend on the type of equation. As initial condition for the numerical solution, use the exact solution during program development, and when the curves coincide in the animation for all times, your implementation works, and you can then switch to a constant initial condition: $$u(x,0)=T_{0}$$. 0 & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & 0 \\ Despite the fact that the Crank-Nicolson method, or the θ rule with $$\theta=1/2$$, is theoretically more accurate than the Backward Euler and Forward Euler schemes, it may exhibit non-physical oscillations as in the present example if the solution is very steep. # Return the final array divided by the grid spacing **2. π 2 ∂ u ∂ t = ∂ 2 u ∂ x 2. In particular, we may use the Forward Euler method as implemented in the general function ode_FE from Sect. As we’ll see in the next chapter in the process of solving some partial differential equations we will run into boundary value problems that will need to be solved as well. When the temperature rises at the surface, heat is propagated into the ground, and the coefficient β in the diffusion equation determines how fast this propagation is. Boundary value problems arise in several branches of physics as any physical differential equation will have them. In the present case, it means that we must do something with the spatial derivative $$\partial^{2}/\partial x^{2}$$ in order to reduce the partial differential equation to ordinary differential equations. Show that if $$\Delta t\rightarrow\infty$$ in (5.16)–(5.18), it leads to the same equations as in a). … There is also diffusion of atoms in a solid, for instance, and diffusion of ink in a glass of water. The CPROP based approach is extended to a constrained integration (CINT) method for solving initial boundary value partial differential equations (PDEs). 107.170.194.178, We shall focus on one of the most widely encountered partial differential equations: the diffusion equation, which in one dimension looks like, $$\frac{\partial u}{\partial t}=\beta\frac{\partial^{2}u}{\partial x^{2}}+g\thinspace.$$, $$\frac{\partial u}{\partial t}=\beta\nabla^{2}u+g\thinspace.$$. Open Access This chapter is distributed under the terms of the Creative Commons Attribution‐NonCommercial 4.0 International License (http://creativecommons.org/licenses/by-nc/4.0/), which permits any noncommercial use, duplication, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, a link is provided to the Creative Commons license and any changes made are indicated. It reads: The next equation - around grid node 2 - reads: For this one, there is nothing to change. The focuses are the stability and convergence theory. However, these authors prefer to have an ODE for every point value u i , $$i=0,\ldots,N$$, which requires formulating the known boundary at x = 0 as an ODE. T_{j+1}\\ In addition, the diffusion equation needs one boundary condition at each point of the boundary $$\partial\Omega$$ of Ω. Matlab/Octave contains general-purpose ODE software such as the ode45 routine that we may apply. Identify the linear system to be solved. Matrix and modified wavenumber stability analysis, 4. \vdots \\ To avoid oscillations in the solutions when using the RKFehlberg method, the rtol and atol parameters to RKFFehlberg must be set no larger than 0.001 and 0.0001, respectively. Report what you see. The main advantage of this scheme is that it is unconditionally stable and explicit. \vdots \\ \frac{-\frac23 T_1 + \frac 23 T_2}{\Delta x^2} = b_1 + \frac{4}{3 \Delta x}.\], # right-hand side vector at the grid points, Constructs the centered second-order accurate second-order derivative for, matrix to compute the centered second-order accurate first-order deri-, vative with Dirichlet boundary conditions on both side of the interval. # Perform the matrix multiplication of the inverse with the right-hand side. Consistency and monotonicity of the scheme are discussed. A major problem with the stability criterion (5.15) is that the time step becomes very small if $$\Delta x$$ is small. You can print out solver_RKF.t_all to see all the time steps used by the RKFehlberg solver (if solver is the RKFehlberg object). In fact, a large part of the solution process there will be in dealing with the solution to the BVP. T_j \\ T_{nx-2} # Manually set the boundary values in the temperature array. b_{j+1}\\ In mathematics, in the field of differential equations, a boundary value problem is a differential equation together with a set of additional constraints, called the boundary conditions. Solve this heat propagation problem numerically for some days and animate the temperature. Over 10 million scientific documents at your fingertips. With N = 4 we reproduce the linear solution exactly. We have seen how easy it is to apply sophisticated methods for ODEs to this PDE example. for solving partial differential equations. 1 & -2 & 1 & 0 & 0 & \dots & 0 & 0 & 0 & 0 \\ Filename: symmetric_gaussian_diffusion.m. Occasionally in this book, we show how to speed up code by replacing loops over arrays by vectorized expressions. Non-homogeneous Dirichlet boundary conditions. Here, we will limit our attention to moderately sized matrices and rely on a scipy routine for matrix inversion - inv (available in the linalg submodule). Assume that the rod is 50 cm long and made of aluminum alloy 6082. For the left boundary node we need something different. *y(1: end)). Here is the Python code for the right-hand side of the ODE system (rhs) and the K matrix (K) as well as statements for initializing and running the Odespy solver BackwardEuler (in the file rod_BE.py ): The file rod_BE.py has all the details and shows a movie of the solution. What takes time, is the visualization on the screen, but for that purpose one can visualize only a subset of the time steps. Partial Diﬀerential Equations Igor Yanovsky, 2005 2 Disclaimer: This handbook is intended to assist graduate students with qualifying examination preparation. The Odespy solvers expect dense square matrices as input, here with $$(N+1)\times(N+1)$$ elements. For example, halving $$\Delta x$$ requires four times as many time steps and eight times the work. In other words, with aid of the finite difference approximation (5.6), we have reduced the single partial differential equation to a system of ODEs, which we know how to solve. By B. Knaepen & Y. Velizhanina The solution of the equation is not unique unless we also prescribe initial and boundary conditions. 1st order PDE with a single boundary condition (BC) that does not depend on the independent variables Linear PDE on bounded domains with homogeneous boundary conditions # We use d2_mat_dirichlet() to create the skeleton of our matrix. Note how the matrix has dimensions (nx-2)*(nx-2). As we are using a second-order accurate finite difference for the operator $$\displaystyle\frac{d^2 }{dx^2}$$, we also want a second-order accurate finite difference for $$\displaystyle\frac{d }{dx}$$. By using this website, you agree to our Cookie Policy. \vdots \\ The ode_FE function needs a specification of the right-hand side of the ODE system. The initial and boundary conditions are extremely important. The equation is written as a system of two first-order ordinary differential equations (ODEs). However, since we have reduced the problem to one dimension, we do not need this physical boundary condition in our mathematical model. \\ The matrix $$\tilde A_{ij}$$ on the left-hand side has dimensions $$(nx-2)\times(nx-2)$$. Note that all other values or combinations of values for inhomogeneous Dirichlet boundary conditions are treated in the same way. You can then compare the number of time steps with what is required by the other methods. Diffusion processes are of particular relevance at the microscopic level in biology, e.g., diffusive transport of certain ion types in a cell caused by molecular collisions. Around the other grid nodes, there are no further modifications (except around grid node $$nx-2$$ where we impose the non-homogeneous condition $$T(0)=1$$). 0 & 1 & -2 & 1 & 0 & \dots & 0 & 0 & 0 & 0 \\ We remark that a separate ODE for the (known) boundary condition $$u_{0}=s(t)$$ is not strictly needed. The physical significance of u depends on what type of process that is described by the diffusion equation. # We only need the values of b at the interior nodes. Partial differential equations. Introduction to Differential Equation Solving with DSolve The Mathematica function DSolve finds symbolic solutions to differential equations. Our setting of parameters required finding three physical properties of a certain material. Then, taking into account $$T_{nx-1}=0$$, the equation around grid node $$nx-2$$ becomes: If we collect these $$nx-2$$ equations back into a matrix system, we get: The above system is completely closed in terms of the actual unknowns $$T_1,\dots,T_{nx-2}$$. To summarize, the partial differential equation with initial and boundary conditions reads, $$\frac{\partial u(x,t)}{\partial t} =\beta\frac{\partial^{2}u(x,t)}{\partial x^{2}}+g(x,t), x\in\left(0,L\right), t\in(0,T],$$, $$\frac{\partial}{\partial x}u(L,t) =0, t\in(0,T],$$, $$u(x,0) =I(x), x\in\left[0,L\right]\thinspace.$$, x_{0}=0