logo
  • userLoginStatus

Welcome

Our website is made possible by displaying online advertisements to our visitors.
Please disable your ad blocker to continue.

Current View

Biomedical Engineering - Mathematical and Numerical Methods in Engineering

Completed notes of the course - Numerical Methods

Complete course

LEZIONE 1 – 13.09 GENERAL CONCEPTS Why do we study this c ourse for PDE. The final objective of this curse is to give you a feeling of what the MM based on PDE are helpful for the description of phenomena that you will study. All fields of engineering are about PDE. This course will help you in how to approach at new model, and how to use the computer for your help. You have to understand that with mathematical tools you can represent all of your interest. This is the purpose of this course. We have more to ols to analyze hard problems. All the models of continuum mechanics are based on 2 fundamental ingredients: the conservation principles (of mass, momentum, energy, charge…), which are theoretically derived, that you have to translate in mathematical model s, and the constitutive laws where you reconnect, experimentally determined (diffusion, friction, chemical kinetics…) and specific of all various fields. This is called phenomenological approach, because you observe the phenomena and try to learn his laws through this approach. Just a few examples. 1. SOLID MECHANICS: ELASTICITY You just take mass and momentum conservation into 3 dimensions for an object, and you combine it with laws about deformation. The stress is proportional to the deformation , so, if you pull the stress increases like a coil and writing this in a form of 3 dimensional setting (infinitesimal calculus), you obtain the information of solid in 3D, and this is fundamental for the design. 2. FLUID MECHANICS: NAVIER -STOKES is the principle of conse rvation of mass for a fluid that is called computational fluid dynamics (CFD). 3. FLUID STRUCTURE INTERACTION In this picture you can put all things together: the interaction of a solid with a fluid. The combination of mechanics is being the dance under the action of a fluid that generate some forces on the beam. Fundamental biological phenomena are subject to this type of interaction, for example cell in general or bacteria. AN OVERVIEW OF MODELS BASED ON PDEs These equations are the equation of heat balance. U is the temperature  scalar field Φ is the heat flux  vector field The equation 1 describes the variations of temperature in a room (time/space). The temperature varies in this room provided that there is a heat exchange. So, this equation is basically about heat balance or heat exchange. There is a heat going in and out the room: if I open the room there is some cool air and the temperature will decrease. The temperature increases in time if the heat goes in, and in the same way the temperature decreases in time if the heat goes out. This is the balance principle. So, let’s translate it into an equation. Temperature increase means that we are speaking about time by time derivative . Temperature derivative: ������������������= HEAT GOES IN . This is the preliminary mathematical setting of heat balance principle. The only thing that we have to quantify is what does it mean that the heat goes in. what does it mean in mathematical language? Imagine the situation in which we have the heat that g oing in in all the possible directions, in this case we have to introduce the theme of the divergence of the heat f lux (vector field). Divergence of heat flux : ������⋅������ = ��������������+ ��������������+ ��������������, ������ = (�������,�������,�������), ������ = (�������,�������,�������) DIVERGENCE THEORE M The divergence is a measure of how much source there is in each infinitesimal part. ������������������= HEAT GOES IN = −������⋅������. Negative because goes in. ������������������ + ������⋅������ = 0  THE HEAT BALANCE EQUATION. We don’t use the heat balance equation, but the constitutive laws, and by this of heat transfer, called the Fourier Law: ������ = −�������������������, D = coefficient of conductivity  1D ������ = −������������������  2D,3D Now just put things together: ������������������− ������(������������������ )= 0  2D,3D , ������������������− �������(�������������������)= 0  1D => ������������������− �������������� ������= 0 If the temperature is not uniform, it will change according to this type of equation, called the heat equation. PROTOTYPE OF LINEAR PDE With this approach we study different types of PDE, all summarized i n this general equation. The diffusion term is manifested when you put a drop of hot water in the cold water, and you see that it spread off, with the equation of diffusion/viscous term. But also, the heat in a room could be spread off with the conditionate air, with the diffusion equation but also with the advection term. So, we have a combination of all these terms. Another example is the Bacteria model, that study how the bacteria consume the nutrients, like the glucose. This model is based on the diffusion reaction equation. DIFFUSION Image that this picture is the heat in a room, with high and low region. This, according with the equation that we have just formulated, will melt down by diffusion , but also by advection if it moves. It models processes that diffuse without preferential directions of propagation . ADVECTION It models processes that propagates in the direction of a . The shape remains unchanged but moves. The combination of diffusion and advection is a prototype of a very important class of phenomena that are dynamics. All fluids are the combination of these two. DIFFUSION -REACTION If we combine diffusion and reaction, we obtain that is called pattern for mation, equation that describes all patterns. It models local reaction processes . For example, studying this equation we can interpret the pattern of the skin of the animals. The shape of the pattern will depend on the size of the animal, because the equat ion acts on a smaller or bigger size giving rise a different types of solutions, and the pattern is the visualization of the solution of this equation. For example, snakes can’t have dots because the expect ratio is too long and there is not enough room on the snake to develop oscillations in circular direction. So, snakes could have stripes, but not dots. The goat is an example of the visualization of cosine function, because half of it is black (when cosine is positive), and the other one is white (becaus e cosine is negative). That’s very interesting… but is it useful? In the reality, the same principle determines the engineering devices here. In these pictures there are stents born from the mathematical models. They are metal devices that were implanted percutaneously in arteries after angioplasty for dilate lumen in case of arteriosclerosis. They have to be flexi ble to navigate to the point , also quite stiff to support the arteria wall and have to be less invasive to not interact too much with the blood flow. Sometimes these devices are in polymer and able to release some drugs, chemical agent that helped the arte ry to follow the right evolution. They have physical models that we have talk about: mechanic for the stent analysis, fluid mechanics for the interaction with the blood flow and mass conservation for the drug release. MAIN TOPICS There are two parties: th e numeric in blue and the analysis in red; the interaction of two is the whole thing that you have to understand. Finally, in practice we have the numerical labs , in which with Matlab we see the simple problem and the behavior of the equation of numerical method. Review of some topics of numerical analysis . Solid mechanics problem : the discretization of linear PDE (for example diffusion or the heat equation or the elasticity equation) ends up with definition and finally the solution of large possibly sparse linear system of equation . Find vector x:Ax=b, where the dimension of the matrix is n and is very big, as 10^6, while x ∈ ℝ������, A ∈ℝ������������� . Most of the computational time is spent in the solution of the linear system , because most of the computationa l time is spent in process of resolving linear system. NUMERICAL SOLUTION OF LINEAR SYSTEM The black line is increased of computing power  MOORE’S Law. From this curve we can see that there was an efficient improvement since 1970. The red one is the effic iency of numerical algorithms , that is independent form the computational power but has a comparable rate. Th ese improve with passing of time, until reach very big problems , like Earth simulator for understanding the trends for climate change. So over the years the numerical algorithms became more difficult . Idea (Gauss): Gauss elimination method . Is a direct method for solving linear system: in a finite number of operations we can find the solution of any well posed linear system. The exact so lution is x. If I take a system Ax = b of dimension n, how many operations shall we expect to do before the ending up the algorithm? we expect the exact solution in n^3  ������(�3). E.g.: �= 10 6, ������(�3)= 10 18. The method is a s equence of operations that sta rting from a simpl e linear system, in this case 3x3 system, tries to transform the system in a triangular system with the same solution . The advantage is that we can apply t he backward substitution , and this only cost n^2 operations . Backward substitut ion We can interpret it in an elegant form. Doing GEM is equivalent to calculate two triangular matri ces L,U (lower and upper triangular) . Such that they provide the factorization of matrix A; a factorization is multiplicative decomposition of the matrix  A = LU. Why this decomposition is so helpful? LU is a new big triangular, so is very helpful to solve Ax = b LU DECOMPOSITION: We can transform the original system into the combination of two triangular system . Y is the auxiliary unknow, so Ux = y; so, for the first we use a forward substitution and for the second a backward substitution. At this point, we have an interesting conceptua l step. All the method that followed here are iterative methods (after 1970) . The idea of iterat ive method is to approximate the solution of A x = b using a sequence of vectors  ��∈ ℝ������; such that they satisfy the fundamental property of convergence. For computational needs stop when: k>1 , typically n = 10^6  K(A) = 10^12. In this case the contraction factor of this formula is very close to 1, meaning that t he convergence is very slow. So, the solution could be preconditioned. PRECONDITIONING : P is the preconditioner. ������−1������� = ������−1�, then �(������−1������) �(��+1)= �(��)+ ������� �′(��)+ 1 2�������2�′′(��)  TRUNCATED TAYLO R EXPANSION (to the 2 nd order) We are certain that exist a point such that this formula is correct: ∃��∈(��,��+1). We can replace the formula above, and we obtain: �′(��)− 1 ������� [(�(��)+ �������� ′(��)+ 1 2�������2�′′(��))− �(�������)] = ?  − 1 2������� �′′(��)  error is linear ly proportional to ������� .  first order scheme. Recap : If you do same things to B and C scheme, you obtain this: B and F schemes are equal (1 st order), while C scheme is a 2 nd order. APPROXIMATION OF SECOND ORDER DERIVATIVES Why is it important? For example, the heat equation in one dimension is: ���− ��������� �= 0, D = diffusion operator is relate d to the secon d derivative. In mechanical problems, more precisely to the elastic deformations in 1D: −���� �= �, u is a displacement, f are forces. The question is: given a function y(t) approximate y’’(t) ; the idea is: ⅆ ⅆ�(ⅆ ⅆ��(�)), approximate the outer derivation (fuori parentesi) with a center scheme based on a half -step. It means that: ⅆ ⅆ��′(�������)= �′(�ⅈ+12)−�′(�ⅈ−12) ������� .  ⅆ ⅆ�(ⅆ ⅆ��(�)) = ⅆ ⅆ��′(��). The second idea is using a center scheme for the inner derivatives, based on a half -step: Y’(��+12)≈ �(�ⅈ+1)−�(�ⅈ) ������� ; Y’ (��−12)≈ �(�ⅈ)−�(�ⅈ−1) ������� . In conclusion: �′′(��)≈ 1 ������� (�(�ⅈ+1)−�(�ⅈ) ������� − �(��)−�(�ⅈ−1) ������� )= �(�ⅈ+1)−2�(�ⅈ)+�(�ⅈ−1) �������2 . Recap: derivation of the for mula. Back to differential equation (First Order Cauchy Prob lem): �′(�)= �(�,�(�)). A numerical scheme is a new problem that defines the numerical solution, meaning that: ��≈ �(�������), the forward numerical scheme at time ti, �ⅈ+1−�ⅈ ������� = �(��,��). This formula defines the sequence or the discrete function ui, and this is a possible candidate for numerical scheme. The backward scheme at ��+1: �ⅈ+1−�ⅈ ������� = �(��+1,��+1)  backward approximation in time with respect to ��+1. This is an approximation of a mathematical problem with another problem. EU LER SCHEMES We have to study the main concept that allows us to characterize the three schemes. This comparison is like hav ing different recipes (F,B,C), which one is ‘good’, which one is the best. We would know what the best scheme is to use. ERROR ANALYSIS: we have to define what is the error of each numerical scheme. The error is the distance/difference between the exact solution and the numerical solution. For example, for the Cauchy problem, the exact solution was named y(t), the numerical solution un. So , the error at time tn is: �������= |�(�������)− �������|, ∀������= 0,1,2,⋯ ������. This is a very important concept about a numerical scheme. The exact value of the error is UNKNOW, meaning that otherwise �(�������)= �������+ �������. This is not true in general, because usually is unknow . What is important is this: The information about how the error scales with respect to the discretization parameter, more precisely if the error goes to zero asymptotically with the discretization parameter, which in this case is ������� . In mathematical te rms we want to determine if �������= ������(�������������)  big O (o grande), ������(������)= ������⋅������, where C is independent of a. Here we want to understand if the error is proportional to ������� to the power P, where P is a non -negative index (to have convergence P>0), and P is named order of scheme. Every time we consider a numerical scheme, we want to know the behavior of the error. Examples of the Cauchy problem discretized by Euler method the error can be combined by two factors: - The intrinsic error of one single step of the scheme - The combination (the same) of all previous error The most important error is the intrinsic of one step. How do we define it? LOCAL TRUNCATION ERROR (LTE) Imagine that at time n, the numerical approximation is coincided with the exact solution.  hy pothesis: �������= �(�������). Then we take one step of the scheme. Thanks to the dashed line, we obtain an approximation at time �������+1, which is indeed not the exact solution, because by taking one step of the scheme we didn’t do an exact operation with respect to the actual Cauchy problem. So, we have in this case a special form of the numerical approximation, called �������+1 ∗ , which is the numerical approximation obtained when the approximation of previous step is exact: �������= �(�������). NO ACCUMULATED ERROR. The int rinsic error of one step is just the difference between the analytical solution at time �������+1 and this special numerical solution that we obtain when we cancel the previous accumulated error: �(�������+1)− �������+1 ∗ . Using the FE in the canonical form �������+1= �������+ �������� (�������,�������): �������+1 ∗ = �(�������)+ �������� (�������,�(�������)). The intrinsic error of one step is �(�������+1)− �(�������)− �������� (�������,�(�������))  ������� ⋅������������������ = �(�������+1)− �(�������)− �������� (�������,�(�������)) This is for FE. The main concept is that the error can be decomposed in 2 parts: the intrinsic error of one step and the error that has been accumulated from the beginning. Exercise Find the order of the LTE of the FE scheme with respect to ������� . Let’s start from the expression of LTE. ������� ⋅������������������ = �(�������+1)− �(�������)− �������� (�������,�(�������)) = ������(�������������). We take the truncated Taylor expansion: �(�������+1)= �(�������)+ ������� �′(�������)+ 1 2�������2�′′(�������)  �������∈(�������,�������+1). Substituting we obtain: ������� ⋅������������������ = ������� �′(�������)+ 1 2�������2�′′(�������)− �������� (�������,�(�������)). Y(t) : y’(t) = f(t,y(t))  the solution of the Cauchy p roblem  �′(�������)= �(�������,�(�������)) So, ������� ⋅������������������ = ������� �′(�������)+ 1 2�������2�′′(�������)− �������� (�������,�(�������))  ������� ⋅������������������ = 1 2�������2�′′(�������). This means that ������������������������ = 1 2������� �′′(�������)= ������(������� ). Recap: for the error analysis of the finite difference scheme, we will first of all analyze the LTE of the scheme and then study how the LTE at each time step combines to give the total error, ������� . CLASSIFICATION Some theorems help us to do a classification between one scheme and the others, making a distinction, and to understand what scheme is the best. Very important types of criteria that we will use in general to make distinction between the schemes are the follow:  Implicit/Explicit scheme s: Forward Euler is explicit (we can calculate the numerical solution at new time as correction of the numerical solution at the previous time – conditionally stable ). Backward Euler is implicit (where �������+1 can be expressed as a correction of the previou s calculated numerical solution, but the increment depends on the new solution at new time – nonlinear problem – unconditionally stable).  One step /multistep: for PDEs mostly see the one step scheme  Order of the scheme: what we discuss previously for LTE o r GTE. I s the parameter P such that the error is �������= ������(�������������)  ������������������������ = ������(�������������). This is very important because the higher P is, better the scheme is. The question now is: how do we analyze the total error, �������, which is a combination of intrinsic error at every step ? In general in all numerical scheme, but in particular the general framework for finite difference scheme is the so called EQUIVALENCE PRINCIPLE (or the LUX equivalence principle). EQUIVALENCE PRINCI PLE The main principles that define a numerical scheme are these: The final objective is the convergence of the error: �������→ 0 as ������� → 0. The convergence is a complex property that depends on one hand on the intrinsic error and on the other hand on the accumulation of the intrinsic errors. So, the equivalence principle explains that these 2 effects are controlled by 2 different simple properties, the consistency and the stability. Consistency: related to the intrinsic error of the scheme Stability: relate d to the accumulation of error along the time step Basically, the concept of the equivalent principle is quite simple and is that if the intrinsic error of the scheme is small and proportional to ������� , and if along the sequence of time step, the error of e ach time step doesn’t accumulate too much, I can co mbine these 2 elementary properties in order to infer the convergence of the scheme. So, the fundamental importance of this principle is that if I want to study the convergence, I have to first determine if the consistency of the scheme is true and if the stability of the scheme is true. Then, by the general framework of the principle we c an derive the convergence. This is the center of all the course, of all numerical analysis. Why stability affects the accumulation of errors? Stability in general is the property such that small perturbations (or small errors) on the initial condition propagate into small errors on all solution, as long the solution evolve s with time . So, small perturbations are the initial conditions and propagate into small errors along the evolution of the system. In the example we have our dif ferential equation: �′= �(�,�), �(0)= �0. We see that the evolution of the system unperturbed is this black line, is the trajectory of the system for the exact value of initial condition. Imagine that now we introduce a little perturbation of the initial c ondition: the trajectory of the system will be affected, but for a stable system the perturbated trajectory will always be close to the original one (red line). This basically explain directly why stability affect the accumulation errors: because errors ar e small perturbation introduced at every time step. The effect of all perturbation will remain bounded, and, in this way, the total accumulated error will remain on the same magnitude of each single intrinsic error. And for this reason, the total accumulat ed error will be also proportional to ������� . So, the stability really guarantees that the errors don’t grow to much along the evolution of the system, and this is the fundamental role of the stability in every numerical scheme. Stability can be split in 2 different more precise property: - Zero -stability : measure the stability with respect the parameter ������� or h when h -> 0 ; a scheme is zero -stable if it is stable for any h -> 0 for a fixed time - Absolute stability: measure the parameter n , the number of time step, when it goes to infinity; a scheme is absolutely stable if it is stable for N,T -> infinity, fixed h For PDE we will talk about stability in general, without this distinction. A consequence of the equivalent principle is that stability combined with consistency, which is about ������������ ������������→ 0 as ������� → 0, imply the fundamental property of convergence: Initial perturbation METODI ANALITICI E NUMERICI PER LE EDP SCALAR AND LINEAR CONSERVATION L AWS – TRANSPORT EQUATION A transport equation, for example , g iven a parameter a (a>0, but not always), given the initial state �0(�), we want to find a function �(�,�) such that ��������(�,�)+ �������������������(�,�)= 0. It means that this function has a bump in a point, and then extend to zero all time. Is equivalent to initial condition for ODE, but the theorem is not just a condition in one point, but in all real line (x – axis). The physical meaning of transport problem is to transport this function almost unperturbed in one direction, in this case to the left because a>0, al ong the space. So, the solution of the problem is the initial bump, that translates along the space with a given constant velocity (a). This is just the effect of pure translation. Slightly more complex models can describe important phenomena such as shock waves in gasses  Burgers equation: ��������+ ��������������= 0. The problem, when it is scalar and linear, is just a pure transport problem. We can combine 2 important effects: transport and diffusion. This is typically an important type of phenomena in many examples, for example fluid mechanics. APPLICATIONS 1. Traffic flow: can be modeled describing the density of cars ; the maximum speed is a>0. We can use the transport equation to describe how distribution of cars moves along the highway, and under a perturb condition we can see a pack of cars that only transl ate. We haven’t considered other aspects, like acceleration or deceleration of cars, that might make the system more difficult. So, shock waves can correspond to the formation of traffic jams. 2. Gasdynamics: how speed aerodynamics. We take a one -dimensional tube with a septum that separate 2 parts of the tube, one with high pressure and the other with low pressure. Then we break the septum, and we see the gas expand very suddenly. 3. Particle transport: for example, the propagation of the spray in the airways t o study how they are transported by the breath into the lungs. These include maybe transport and diffusion. THE FINITE DIFFERENCE METHOD Now, let’s discuss about the derivation of the main schemes for the transport problem. The main concept is to appr oximate the derivatives of the equation ��������(�,�)+ �������������������(�,�)= 0, using discr ete derivatives. The difference between this problem and the Cauchy problem is that they are initial value problem both, but the Cauchy problem has only one independent variabl e (t), and here we have two independent variables (t,x), and the derivative are partial derivatives, with respect time only or space only. In order to introduce the numerical discretization based on final differences, we have to introduce the main concept s for obtaining a numerical scheme. - the computational grid/mesh: a discretization of space and time in two points - the numerical approximation: some formulas for describe derivatives , for example to approximate the derivative in time and the derivative in space - the numerical scheme: to calculate the numerical solution, which is an approximation of the exact solution that we may not be able to calculate analytically NUMERICAL DISCRETIZATION In order to define a computational grid, we place nodes or points along x and along t. h is the space discretization step, sometimes indicated with ������� . τ or ������� is the time discretization step. The objective of discretization is to approximate the equation at the points of the computational grid, defined for example by (�������,�������). ������������� is the numerical solution in (�������,�������).  ������������� ≈ �(�������,�������). The numerical scheme is the algorithm that give s us �������������. APPROXIMATION IN TIME: BACKWARD AND FORWARD SCHEMES Use discrete derivatives on the computational grid to define a scheme . For example, taking the time derivative in one point: �������(�������,�������)= �(�������+1,������������)−�(�������,������������) ������� + ������(������� )  FORWARD IN TIME ��������(�������+1,�������)= �(�������+1,������������)−�(�������,������������) ������� + ������(������� )  BAC KWARD IN TIME APPROXIMATION IN SPACE: CENTERED AND ONE -SIDE SCHEMES �������������(�������,�������)= �(�������,������������+1)−�(�������,������������−1) 2������������ + ������(�������2) CENTERED IN SPACE �������������(�������,�������)= �(�������,������������)−�(�������,������������−1) ������������ + ������(������� )  ONE -SID E IN SPACE We can put together to obtain these: Backward in time, one -sided is not often used, because is not needed. For the criteria of stability, we will see that forward in time, centered is unstable, or conditionally stable under s evere restriction; forward in time, one -side is conditionally stable under moderate condition; backward in time, centered is unconditionally stable. (1)  ������������ ������������������= ������(������� ,�������2); (2)  ������������ ������������������= ������(������� ,������� ); (3)  ������������ ������������������= ������(������� ,�������2); To analyze the differences between these schemes, we have to talk about the stencils STENCILS Is the collection of nodes that are involved with the computational of the scheme . For example, if you take in blue the forward, centered what we obtain that the stencil is com posed by the points in the picture. If you consider the forward, one -sided (in green). While, for backward, centered (in red). The forward in time schemes are explicit schemes, because given the solution in the previous time step, we can calculate the so lution at the new time step by simple operations. While, for the backward scheme, all the points at the time �������+1, are coupled together and consist of numerical unknows and this makes it an implicit scheme. Replace ������������� to u (�������,�������) at any node The objective of today is had another few examples and cast all of them in general framework, and finally move towards the analysis of the scheme to understand their properties and the convergence. Computational standpoint: how to translate / implement the scheme on a computer and solve it? The behavior of the Implicit and Explicit Euler is easy to study if we l ook at the ma trix /vector formulation of the scheme. It means that we represent the unknown numerical solution as a vector, which component is one of the values of numerical solution at the node ������������, and we represent the schemes using a matrix vector multiplication. This is important for the implementation of the scheme, especially in matlab. For example, we start with FE/C: �������������+1= �������������− ������������� 2������������ (�������+1������ − �������−1������ ). Now we want to collect in a single vector all the unknows that represents the approximations of values at the nodal points. �������+1 = | �1������+1 �2������+1 �������������+1 |  approximation of the nodes. In a similar way we can define �������. This is a formula that connects ������������� to �������, a linear combination of components of �������. In the end we want to express the numerical scheme as this expression: �������+1 = �ℎ������*�������, τ = ������� and h = ������������ . �ℎ������is a matrix. For each scheme determine �ℎ������. Let’s see how we do it.  example about the upwind scheme. Vector/matrix form of an explicit scheme FE/C: �������������+1= �������������− ������������ 2ℎ (�������+1������ − �������−1������ ).  �������+1 = ������� −������������2ℎ�*�������, B is the matrix. B in the central diagonal has 0, because there is not i coefficient here. In the lower diagonal we put the coefficient i -1, and there is -1 because of the negative sign . While, for i+1, there is +1 because of the two negative sign in the last formula. Diagonal that corresponds to the index i Lower diagonal Upper diagonal So, �������+1 = ������� −������������2ℎ�*������� = (������− ������������ 2ℎ �)∗ �������. We can conclude saying that, with the final notation, for FE/C: �ℎ������ = (������− ������������ 2ℎ �). This is the typical structure of an explicit scheme: �������+1 = �ℎ������*�������. Vector/matrix form of an implicit scheme Example: backward Euler/centered , BE/C �������������+1= �������������− ������������ 2ℎ (�������+1������+1− �������−1������+1)  �������������+1+ ������������ 2ℎ (�������+1������+1− �������−1������+1)= �������������. �������+1 + ������������ 2ℎ �*�������+1 = �������.  the general form is this: �������+1 = �ℎ������*�������  (������+ ������������ 2ℎ �)∗ �������+1 = �������; �ℎ������ =  (������+ ������������ ℎ�)−1  the inverse This expression means that the BE/C involves the solution of a linear system at each time step, and this indeed the typical behavior of the implicit scheme. FINITE DIFFERENCE SCHEMES Let’s summarize this introduction. We have introduced the first 3 schemes for linear conservation law that can be classified in different ways: forward in time (explicit behavior) or b ackward in time (implicit behavior), one -sided or centered (relative to space approximations). We can add some new options to these. To understand how, let’s start from the upwind scheme. Upwind Scheme Why is it called upwind? Imagine that I want to appro ximate a transport problem: a > 0  transport parameter positive, is the velocity of transportation. In the case positive, the translation is from the left to the right. For example, we translate the initial state of this form, with speed a, in the direc tion of axis x to find the true solution of the problem. So, the meaning of the upwind scheme is to break the symmetry of the center of the approximation, in a way such that the numerical scheme reproduces the direction of transport. It means that if I hav e to approximate the solution in the blue point in the first graphic, and I want to break the symmetry don’t using the two nodes at ������������−1 and ������������+1, I decide to mimic the directional transport, so I will use the point at same location of the previous ti me, and the point that is upwind to xi. It means that, by coupling these 2 points, I basically say that the numerical information in the scheme travels in the same direction as the physical information that travels with positive speed a. So, the numerical information travels in the same direction of a. let’s see what happen in the case of negative a: a 0: �������������−�������−1������ ������������ - a < 0: �������+1������ −������������� ������������ in the end, the UW scheme has a multiple form: (�������������− �������−1������ ) �������������+1= �������������− ������������ ℎ  the general form of UW. We know why it is called upwind. (�������+1������ − �������������) Sometimes, is not very practical to write the scheme in this form, so there are two alternative formulas. Can we write this scheme with one single formula? (Combining the two options) The answer is yes, we can: �������������+1= �������������− 1 2������������ℎ(�������+1������ − �������−1������ ) + 1 2|������|������ ℎ(�������+1������ − 2�������������+ �������−1������ ). Are these two formulas equivalent? The answer is yes. The following is the verification. - Let’s take a>0: |������|= ������. After the calculation we obtain: �������������+1= �������������− ������������ℎ(�������������− �������−1������ ) - Let’s take a |�������������|. The coeff icient is the value of p. If p increases , we wait to more the larger nodal values of the solution. What happen if we take p to ꚙ? This is the infinitive norm: ‖������������‖������,∞ = h max������ |�������������|. Stability of FE/C,BE/C,UW,LF,LW We will prove 3 fundamental properties of the numerical schemes. 1. PROPERTY 1 : is about all the forward schemes, with exception to FE/C. The first property is somehow the standard case for forward schemes. If the CFL condition is satisfied, the UW,LF,LW schemes are strongly stable on the discrete norm with p=1 . For example: ‖������������+1‖ℎ,1≤ ‖������������‖ℎ,1. We must prove this property for UW and LF. Stability of the UW scheme In the case of a>0: �������������+1= �������������− ������������ (�������������− �������−1������ ) (one sided form). ∑ |�������������+1| ������ , the thesis is to prove that is ≤ ∑ |�������������| ������ . The only way to proceed is to use the scheme and s ubstitute on the left -hand side, but before doing that just a little remark. The summation ∑������ goes fro m ⅈ= −∞ ������� ⅈ= +∞. So, I’m not considering the effect of the boundary. This is impossible in practice, in the application of a real scheme, because we always have to cut the domain somewhere on the left and on the right here. Imagine that in this case we are taking the entire real axis domain ( ℝ). Now let’s prove the thesis by substituting on the right -hand side the scheme. We obtain ∑ |�������������+1| ������ = ∑ |�������������− ������������ (�������������− �������−1������ )| ������ = ∑ |(1− ������������ )�������������+ ������������ �������−1������ | ������ . Now it’s the time to apply the CFL condition: |������|������≤ 1, but we are in the particular case where a > 0. So, we have to reformulate the CFL condition: ������������ ≤ 1→ (1− ������������ )≥ 0, ������������ ≥ 0. So, we can now continue our analysis by saying that I can split the absolute value in 2 parts, saying that this is ≤ ∑ |(1− ������������ )�������������| ������ + ∑ |������������ �������−1������ | ������ , and now by the CFL  (1− ������������ )∑ |�������������| ������ + ������������ ∑ |�������−1������ | ������ . Finally, we use the fact that we are taking the infinite summation: ∑ ������� ������ = ∑ �������−1 ������ , so, the shift of 1 index is not relevant. So, I can basically combine these two terms into one, by saying that is ≤ (1− ������������ + ������������ )∑ |�������������| ������ = ∑ |�������������+1| ������ . We proved that ‖������������+1‖ℎ,1≤ ‖������������‖ℎ,1, which is strong stability fo r UW. This is a summary. Stability of the LF scheme The proof goes along the same way as for UW. In this case a > 0 �������������+1= 1 2(�������+1������ + �������−1������ )− 1 2������������ (�������+1������ − �������−1������ ). ℎ∑ |�������������+1| ������ = ℎ∑ |1 2(�������+1������ + �������−1������ )− 1 2������������ (�������+1������ − �������−1������ )| ������ . Then, apply the CFL: ������������ ≤ 1. In this picture, ������ = ℎ. (esercizio per casa  fare da solo lo svolgimento completo) 2. PROPERTY 2 : the BE/C s cheme is unconditionally strongly stable in the discrete norm with p =2. For example: ‖������������+1‖ℎ,2≤ ‖������������‖ℎ,2. 3. PROPERTY 3 : is not very used, because the FE/C scheme is not strongly stable. The scheme is conditionally stable under the condition (������� ≤ (������������ ������)2), more restrictive than CFL. In fact in this case ������� ≤ (������������ ������)2is much restrictive than ������� ≤ (������������ ������), because ������������ is very smaller than a, as in this example: Stability of the BE/C scheme by the energy method �������������+1= �������������− 1 2������������ (�������+1������+1− �������−1������+1), or equivalently [�������������+1− �������������+ 1 2������������ (�������+1������+1− �������−1������+1)]∗�������������+1. This is related to the fact that the property works in the two norm . So, we want to find an expression similar to this: ∑ |�������������+1|2 ������ . Manipulating the previous expression, we obtain this: (�������������+1− �������������)�������������+1+ 1 2������������ (�������+1������+1− �������−1������+1)�������������+1= 0. Now w e apply some algebraic identity to the first term: If we substitute B and A the values of numerical solution, we get (�������������+1)2− (�������������)2+ 1 2(�������������+1− �������������)2. If you continue in this derivation, what do we get for our scheme is the following: (�������������+1)2− (�������������)2+ 1 2(�������������+1− �������������)2+ 1 2������������ (�������+1������+1− �������−1������+1)�������������+1= 0. Take the summation ℎ∑������ : ‖������������+1‖ℎ,22 − ‖������������‖ℎ,22 + ℎ∑ (�������������+1− �������������)2+ ℎ 2������������ ∑ (�������+1������+1− �������−1������+1)�������������+1 ������ ������ . The last term is interesting because is a telescopic summation. It means that, when you do a summation, the terms cancel 2 by 2. So, the terms n+1 cancel with the terms n. So, all summation is equal to zero. And the whole expression as well, equal to zero. So, what we have proved is that, in conclusion: ‖������������+1‖ℎ,22 + ℎ∑ (�������������+1− �������������)2= ‖������������‖ℎ,22 ������ . We also say that the circled additional term is positive. As a result, we obtain what we needed: ‖������������+1‖ℎ,22 ≤ ‖������������‖ℎ,22  strong stability of BE/C. The objectives of today are 2: a) conclude the chapter about scalar conservation laws; b) move on with the new chapter of heat equation The forward Euler scheme, using an approach similar to the backward Euler scheme, is indeed conditionally stable, but the stability condition that characterize it is very restricted. So, for this reason, this scheme is not very used in practice. To use it, we will need to use a time step which is very small, and this isn’t very convenient computationally. In conclusion, we say that we have three different families, according with the criteria that we have introduced. We have forward in time schemes, which are the most used in this case and are explicit schem es. In this category we have UPWIND, LF and LW. UPWIND and LF are first order, while LW is a second order with respect to ������������ and ������� . This last is conditionally stable, because subject to the CFL condition. Then, we have the backward in time schemes: EI C, first order and implicit, because we have to solve the linear system at each time step. This is unconditionally stable. So, this is the typical scheme that would be useful to model a very long time, using a large ������� . Note: we derived and used these s chemes fo r linear problem ��������+ �������������������= 0; �(0,������)= �0(������) , basically the trivial scalar transport equation, which admits an exact solution through the method of characteristic, but the same schemes can be used for general conservation laws, such as: - Burge r equation: ��������+ ��������������= 0; transport equation where the magnitude of velocity is the solution itself  important for fluid dynamics. - Traffic equation: similar to the Burger’s one. THE HEAT EQUATION The heat equation is somehow completely opposite to the transport equation. So, we have two main phenomena that we will be at the basis of the models that we will study. On one hand we have a transport phenomenon, and in the other one we have diffusion phe nomenon. Usually, in nature we find them combined. For example, if we look at the smoke coming out of a chimney, we see that there is a wake, so the smoke is transported by the air or wind, but at the same time the wake become larger and larger due to diff usion. So, in nature we cannot separate one from the other, but in the models we can isolate them. At the beginning we have studied the transport problem, that corresponds to the action of the winds. Now we study the diffusion problem, that basically corre sponds to the pure diffusion and how this transforms a given function, that is the initial state; how diffusion affects, through physical principle the initial state of a system. This is an initial boundary value problem: we have an initial state and at th e boundary of this domain we have a boundary condition , that in mathematical words is called the digital type boundary condition. What happen if we change the boundary condition? (video slide) If we switch to Newman boundary condition, they prevent the co mmunication of the system with the exterior. So, in this case we can assume that the system is a closed system. So, what we are doing is taking one dimensional domain, basically a bar isolated at the two extrema, assigning a given temperature distribution at the initial time, which corresponds to how the temperature equilibrates into the bar. And, because the bar is isolated, the temperature distributes ideally: in other words, the energy of the system distributes equally, uniformly along the system. This i s an example of energy conservation. So, at the beginning the thermal energy is not evenly distributed and equilibrium is perfectly flat. But at any time, you have the same amount of energy in the system. These are all possible solution of the heat equatio n, and the numerical scheme that we will develop are perfectly adequate to model these effects. The problem The heat equation problem is an initial boundary value problem, that can be formulated by stating the heat equation, in this case by two -dimension al notation. So basically, we have the time derivative of function u, that in physics can be related to the temperature, minus the thermal diffusion parameter D, times the Laplace operators appl ied to the temperature field. This is the heat equation, that can be equal to the forcing term, that in this case is a heat source, distributed over the domain. We may have to close the problem to make the problem well post: we have to specify an initial cond ition and possible boundary conditions. We can either have Dirichlet, Newman , or radiation or Robin boundary condition. In one -dimensional problem we can combine these conditions in one extreme, and in the other extreme we can choose any combination. In an y case the problem would be well post. Derivation of the numerical scheme The heat equation is a partial differential equation that depends on time and space: ��������− ������������������ �= � (in this case, the diffusion operator = 1). ������∈(0,1), t > 0. Obviously, this is a time dependent problem, discretized with Finite Differences. And obviously about the computational grid we define some nodes in space: �������= ⅈ������������ = ⅈℎ; and we define also some steps in time: �������= �������� = ������� . Finally, we obtain this comp utational grid. The points where we calculate the numerical solution are only the internal points of the grid. On x0 and xN, the numerical solution is defined by the boundary conditions. The easiest case is Dirichlet boundary condition, that means that the solution is equal to 0 for any time at node x = 0 or x = 1. This implies that �0������ and ������������� are both equal to zero. The numerical approximation consists of calculating the numerical solution ��������= �(�������,�������). This is the main objective of what we want to do at every node of the computational grid. Derivation of the scheme : ��������− ������������������ �= � - Space discretization: ������������������ �(�������,�������)= �(�������,������������+1)−2�(�������,������������)+�(�������,������������−1) ������������2 + ������(ℎ2)  second order scheme Let’s wait for the time discretization. We can formulate an intermediate problem that is called the semi -discrete problem, and is an intermediate step which consist in this: for any node �������, find a function �(�,�������) that is continuous in time and discrete in space, such that satis fy this equation ��������(�,�������)− 1 ������������2(�(�,�������−1)− 2�(�,�������)+ �(�,�������−1))= �(�,�������). (1) Now we introduce the matrix/vector formulation: �(�)= [ �(�,������0) �(�,������1) . . �(�,������������)] , where in our case x0 = 0 and xN = 1. So, I can introduce a vector: �(�)= [ �(�,������0) . . �(�,������������) ] If you use this formulation, the previous equation (1) becomes like this: ⅆ ⅆ��(�)+ 1 ������������2������� (�)= �(�). This is not a partial differential equation, but just a standard differential equation for the vec tor U(t), which corresponds to the first term, and is a combination of values at varying i; this combination of values at varying i can be written multiplying the vector U by a matrix with constant coefficients. The matrix A corresponds exactly to this com bination of coefficients. So, we can also write the expression of A: So, what we immediately see through this approach is that this semi - discrete problem for the heat equation becomes a system of ordinary differential equation with constant coefficient, be cause the matrix A isn’t dependent on U(t). (2) We can use any scheme for the discretization of ODEs (in particular for the Cauchy problem) to discretize problem (2). Problem (2) is called the semi -discrete equation, because we just discretize only the space , but not yet the time. Recap: U(0) is the initial state (Cauchy problem) which is obtained by sampling the initial state �0(������) at the nodes �������. Now, we have to use the numerical schemes for ODEs for the discretization of semi -discrete problem. - Forward Euler - Backward Euler - Crank -Nicholson (trapezoidal) Recap of semi -discrete scheme: - Time discretization of ������̇(�)+ �������������(�)= ������(�); ������(������)= {�������(�������)} What happen if we use a Forward Euler scheme in time ? We look for the solution at time n+1, �������+1, using the solution at time tn, �������. The reference time is �������. ⅆ ⅆ��(�������)= �̇(�������)≈ �(�������+1)−�(�������) ������� . (≈ is an approximation in time). The whole scheme b ecomes: �������+1−������� ������� + ������ℎ�������= �(�������). This is a forward Euler scheme applied to the heat equation, or equivalently the forward Euler centered scheme, center in space forward in time. is the vector that collects all the nodal values of the numerical solution at time n. I don’t include U0 and Un, because of the boundary condition. What happen if we use a Backwar d Euler scheme in time ? The reference time is �������+1. I have to do this: ⅆ ⅆ��(�������+1)= �������+1−������� �������  �������+1−������� ������� + ������ℎ�������+1= �(�������+1). This is the Backward Euler/centered scheme We can now write them in the canonical form, most usual, but multiplying by ������� : FE/C: �������+1= 񘣴