Introduction to Mathematical Modeling and Partial Differential Equations PDF

Introduction to Mathematical Modeling and Partial Differential Equations Mayur J. Bhamborae [email protected] December 2024 What is Mathematical Modeling? Mathematical modeling is the process of using mathematical concepts and language to describe, analyze, and predict real-world phenomena. It involves translating real-world problems into mathematical frameworks that can be analyzed using mathematical tools and techniques. The Modeling Process The mathematical modeling process typically follows these steps: 1. Problem Identification: Clearly define the real-world problem 2. Assumptions: Make simplifying assumptions about the system 3. Formulation: Translate the problem into mathematical terms 4. Analysis: Solve the mathematical problem 5. Interpretation: Translate results back to the real world 6. Validation: Test the model against real data 7. Refinement: Improve the model based on validation Types of Mathematical Models Based on Nature of Variables 1. Deterministic Models: Models where outcomes are precisely determined by known relationships among states and events. dP P = rP (1 − ) dt K (Logistic growth model) 2. Stochastic Models: Models that incorporate random variables and probability distributions. dSt = µSt dt + σSt dWt (Geometric Brownian Motion) Based on Time Dependence 1. Static Models: Describe relationships that don’t change with time. F = ma 2. Dynamic Models: Describe how systems change over time. d2 x + ω2x = 0 dt2 1 Based on Mathematical Structure 1. Linear Models: Systems where output is proportional to input. y = mx + b 2. Nonlinear Models: Systems with more complex relationships. dx = x2 − x3 dt Important Model Types Differential Equation Models Used to describe continuous change: d2 x dx m 2 + c + kx = F (t) dt dt (Mass-spring-damper system) Discrete Models Used for systems that change in steps: xn+1 = rxn (1 − xn ) (Discrete logistic map) Statistical Models Used to analyze data and make predictions: y = β0 + β1 x + ϵ (Linear regression model) Game Theory Models Used to analyze strategic interactions: ui (si , s−i ) (Utility function in game theory) 2 Applications Mathematical models are used in: Physics: Motion, forces, fields Biology: Population dynamics, disease spread Economics: Market behavior, financial systems Engineering: Control systems, signal processing Neuroscience: Neuron modeling, EEG Source Analysis Climate Science: Weather prediction, climate change Social Sciences: Population growth, social networks Limitations and Considerations All models are approximations of reality Models require careful validation The choice of model depends on: – Available data – Required accuracy – Computational resources – Time constraints Differential Equations A differential equation (DE) involves functions of a single independent variable and their deriva- tives. These equations describe how a quantity changes with respect to a single variable, typically time or space. Key Characteristics of Differential Equations Functions depend on one independent variable Use ordinary derivatives (denoted by d dx or d dt ) Describe rate of change with respect to one variable 3 Examples of Differential Equations 1. Newton’s Law of Cooling: dT = −k(T − Ta ) (1) dt where T is temperature, t is time, k is a constant, and Ta is ambient temperature. 2. Simple Harmonic Motion: d2 x 2 + ω2x = 0 (2) dt where x is displacement and ω is angular frequency. Partial Differential Equations A partial differential equation (PDE) involves functions of multiple independent variables and their partial derivatives. These equations describe how a quantity changes with respect to multiple variables simultaneously. Key Characteristics of PDEs Functions depend on multiple independent variables Use partial derivatives (denoted by ∂ ∂x or ∂ ∂t ) Describe rates of change with respect to multiple variables Examples of PDEs 1. Heat Equation: ∂u ∂ 2u =α 2 (3) ∂t ∂x where u is temperature, t is time, x is position, and α is thermal diffusivity. 2. Wave Equation: ∂ 2u 2 2∂ u = c (4) ∂t2 ∂x2 where u is displacement, t is time, x is position, and c is wave speed. Key Differences 1. Number of Variables: DEs: One independent variable PDEs: Multiple independent variables 2. Derivative Notation: DEs: Use dx d notation PDEs: Use ∂x notation ∂ 3. Solution Complexity: DEs: Generally simpler to solve PDEs: Often require more sophisticated solution methods 4 Real-World Applications Differential Equations: – Population growth – Radioactive decay – Simple pendulum motion Partial Differential Equations: – Heat conduction in materials – Wave propagation – Fluid dynamics – Quantum mechanics – EEG Source Analysis The Nabla (∇) Operator The nabla operator, denoted by ∇, is a vector differential operator represented in Cartesian coordinates as: ∂ ∂ ∂ ∇= , , ∂x ∂y ∂z This operator forms the foundation for defining several important vector calculus operations. Gradient The gradient of a scalar field f (x, y, z) is denoted as ∇f or grad(f ). It represents the direction and magnitude of the steepest increase in the scalar field. ∂f ∂f ∂f ∇f = i+ j+ k ∂x ∂y ∂z Properties of Gradient 1. The gradient is perpendicular to level surfaces of f 2. The magnitude of the gradient gives the rate of change in the direction of steepest increase 3. For a scalar field f and constant c: ∇(cf ) = c∇f ∇(f1 + f2 ) = ∇f1 + ∇f2 Example For f (x, y, z) = x2 + 2y 2 + 3z 2 : ∇f = (2x)i + (4y)j + (6z)k 5 Divergence The divergence of a vector field F(x, y, z) = P i + Qj + Rk is a scalar quantity denoted as ∇ · F or div(F): ∂P ∂Q ∂R ∇·F= + + ∂x ∂y ∂z Physical Interpretation Divergence measures the “source” or “sink” strength at a point Positive divergence indicates a source (outward flux) Negative divergence indicates a sink (inward flux) Properties For vector fields F and G and scalar c: ∇ · (cF) = c(∇ · F) ∇ · (F + G) = ∇ · F + ∇ · G Example For F(x, y, z) = (x2 )i + (xy)j + (z 2 )k: ∇ · F = 2x + y + 2z Curl The curl of a vector field F(x, y, z) = P i + Qj + Rk is denoted as ∇ × F or curl(F): i j k ∂ ∂ ∂ ∇×F= ∂x ∂y ∂z P Q R ∂R ∂Q ∂P ∂R ∂Q ∂P ∇×F= − i+ − j+ − k ∂y ∂z ∂z ∂x ∂x ∂y Physical Interpretation Curl measures the rotation or circulation of a vector field The direction of curl is perpendicular to the plane of rotation The magnitude indicates the strength of rotation 6 Properties 1. For vector fields F and G and scalar c: ∇ × (cF) = c(∇ × F) ∇ × (F + G) = ∇ × F + ∇ × G 2. The curl of a gradient is zero: ∇ × (∇f ) = 0 Example For F(x, y, z) = (−y)i + (x)j + (0)k: ∇ × F = (0)i + (0)j + (2)k Important Relationships 1. The divergence of curl is zero: ∇ · (∇ × F) = 0 2. The curl of gradient is zero: ∇ × (∇f ) = 0 3. Double curl identity: ∇ × (∇ × F) = ∇(∇ · F) − ∇2 F The Laplace Operator (∆) he Laplace operator, denoted as ∇2 or ∆, is a differential operator defined as the divergence of the gradient of a function. In Cartesian coordinates for a scalar field f (x, y, z), it is expressed as: ∂ 2f ∂ 2f ∂ 2f ∇2 f = ∆f = ∇ · (∇f ) = + + ∂x2 ∂y 2 ∂z 2 Physical Applications Heat Equation: ∂u = α∇2 u ∂t where u is temperature and α is thermal diffusivity. Wave Equation: ∂ 2u = c2 ∇2 u ∂t2 where c is wave speed. Schrödinger Equation: ℏ2 2 ∂ψ − ∇ ψ + V ψ = iℏ 2m ∂t where ψ is the wave function. 7 Poisson’s Equation: ρ ∇2 ϕ = − ϵ0 where ϕ is electric potential and ρ is charge density. Classification of Partial Differential Equations Classification by Linearity Linear PDEs A PDE is linear if the dependent variable and its derivatives appear only to the first power. General form: ∂ 2u ∂u a 2 +b + cu = f (x, t) ∂x ∂x Examples: Heat equation: ∂u ∂ 2u =α 2 ∂t ∂x Wave equation: ∂ 2u 2 2∂ u = c ∂t2 ∂x2 Nonlinear PDEs Contains the dependent variable or its derivatives with higher powers or within nonlinear functions. Examples: Burgers’ equation: ∂u ∂u ∂ 2u +u =ν 2 ∂t ∂x ∂x Sine-Gordon equation: ∂ 2u ∂ 2u − 2 + sin(u) = 0 ∂t2 ∂x Classification by Order First-Order PDEs Contains only first derivatives of the dependent variable. Example: ∂u ∂u + =0 ∂x ∂y 8 Second-Order PDEs Highest derivatives are of second order. Examples: Laplace equation: ∂ 2u ∂ 2u + =0 ∂x2 ∂y 2 Poisson equation: ∂ 2u ∂ 2u + = f (x, y) ∂x2 ∂y 2 Higher-Order PDEs Contains derivatives of order 3 or higher. Example: Beam equation ∂ 4u ∂ 2u + 2 =0 ∂x4 ∂t Classification by Homogeneity Homogeneous PDEs All terms contain the dependent variable or its derivatives. Examples: ∂ 2u ∂u 2 = ∂x ∂t ∂ 2u ∂ 2u + =0 ∂x2 ∂y 2 Non-homogeneous PDEs Contains terms that don’t involve the dependent variable. Examples: ∂ 2u = sin(x) + t ∂x2 ∂ 2u ∂ 2u + = f (x, y) ∂x2 ∂y 2 Type (Second-Order PDEs) For second-order PDEs in two variables, we consider the general form: ∂ 2u ∂ 2u ∂ 2u A + B + C + lower order terms = 0 ∂x2 ∂x∂y ∂y 2 The classification is based on the discriminant B 2 − 4AC: 9 Hyperbolic: B 2 − 4AC > 0 (commonly describe wave propagation) ∂2u 2 – Example: Wave equation ∂t2 = c2 ∂∂xu2 – Describe wave propagation. Parabolic: B 2 − 4AC = 0 ∂u 2 – Example: Heat equation ∂t = α ∂∂xu2 – Commonly used to model diffusion processes. Elliptic: B 2 − 4AC < 0 ∂2u ∂2u – Example: Laplace equation ∂x2 + ∂y 2 =0 – Describe steady-state or equilibrium problems. – While many familiar elliptic PDEs are time-independent, the classification ”elliptic” itself doesn’t require time-independence. Solving PDEs Analytical Methods Separation of Variables This method assumes the solution can be written as a product of functions, each depending on a single variable: u(x, t) = X(x)T (t) Example for the heat equation: ∂u ∂ 2u =α 2 ∂t ∂x Leading to separated ODEs: T ′ (t) X ′′ (x) = −λα = T (t) X(x) Fourier Transform Method Converts PDEs to ODEs by transforming spatial variables: Z ∞ û(k, t) = u(x, t)e−ikx dx −∞ Particularly useful for linear PDEs with constant coefficients. Laplace Transform Method Transforms the time variable, reducing time-dependent PDEs to spatial ODEs: Z ∞ ū(x, s) = u(x, t)e−st dt 0 10 Green’s Function Method Represents solution as an integral: Z u(x, t) = G(x, t; x′ , t′ )f (x′ )dx′ Where G is the Green’s function satisfying: LuG(x, t; x′ , t′ ) = δ(x − x′ )δ(t − t′ ) Numerical Methods Finite Difference Method Approximates derivatives using discrete differences: ∂u u(x + h, t) − u(x, t) ≈ ∂x h ∂ 2u u(x + h, t) − 2u(x, t) + u(x − h, t) ≈ ∂x2 h2 Finite Element Method Divides domain into elements and uses basis functions: n X u(x) ≈ ui ϕi (x) i=1 Where ϕi (x) are basis functions and ui are coefficients. Spectral Methods Represents solution as a sum of basis functions: N X u(x, t) = an (t)ϕn (x) n=0 Common choices for ϕn (x) include: Fourier series for periodic problems Chebyshev polynomials for non-periodic problems Characteristics Method For first-order PDEs, converts the PDE into ODEs along characteristic curves: dx du = a(x, t, u), = b(x, t, u) dt dt 11 Perturbation Methods For problems with small parameters ϵ: u(x, t) = u0 (x, t) + ϵu1 (x, t) + ϵ2 u2 (x, t) +... Useful for: Weakly nonlinear problems Boundary layer problems Multiple scale analysis Method Selection Choice of method depends on: PDE type (hyperbolic, parabolic, elliptic) Domain geometry Boundary conditions Required accuracy Computational resources The Heat Equation Physical Principles The heat equation is derived from two fundamental physical principles: 1. Conservation of Energy: Energy cannot be created or destroyed 2. Fourier’s Law of Heat Conduction: Heat flux is proportional to the negative temperature gradient One-Dimensional Heat Flow Consider heat flow through a rod of uniform cross-subsectional area A. Fourier’s Law The heat flux q is given by: ∂T q = −k ∂x where: k is the thermal conductivity T is temperature The negative sign indicates heat flows from hot to cold 12 Detailed Derivation Step 1: Consider a Small subsection Take a small subsection of the rod from x to x + ∆x Step 2: Heat Flow Analysis The rate of heat flow through any cross-subsection is: ∂T Q = −kA ∂x Step 3: Net Heat Flow The net rate of heat flow into the subsection is: Net heat flow = Q|x − Q|x+∆x ∂T ∂T = −kA − −kA ∂x x ∂x x+∆x Step 4: Internal Energy Change The rate of change of internal energy in the subsection is: ∂T Rate of energy change = ρcA∆x ∂t where: ρ is density c is specific heat capacity A∆x is the volume of the subsection Step 5: Conservation of Energy By conservation of energy: ∂T ∂T ∂T −kA − −kA = ρcA∆x ∂x x ∂x x+∆x ∂t Step 6: Take the Limit As ∆x → 0: ∂ 2T ∂T kA 2 = ρcA ∂x ∂t 13 Step 7: Final Form Dividing both sides by ρcA: ∂T ∂ 2T =α 2 ∂t ∂x where: k α= ρc is the thermal diffusivity. Physical Interpretation The heat equation shows that: The rate of change of temperature ( ∂T ∂t ) at any point is proportional to the curvature of ∂2T the temperature profile ( ∂x2 ) α determines how quickly heat diffuses through the material The equation is parabolic, reflecting the diffusive nature of heat conduction Physical Setup Consider a vibrating string with: Tension T (assumed constant) Linear mass density ρ (mass per unit length) Displacement u(x, t) from equilibrium position Derivation Steps Step 1: Forces on String Element Consider a small element of string from x to x + ∆x: At point x, tension force is T sin θ1 (vertical component) 14 At point x + ∆x, tension force is −T sin θ2 For small angles: ∂u sin θ ≈ tan θ = ∂x Step 2: Newton’s Second Law For the string element: F = ma ∂ 2u (T sin θ1 − T sin θ2 ) = (ρ∆x) ∂t2 Step 3: Simplification Using small angle approximation: ∂ 2u ∂u ∂u T − = ρ∆x ∂x x ∂x x+∆x ∂t2 Step 4: Take the Limit As ∆x → 0: ∂ 2u ∂ 2u −T = ρ ∂x2 ∂t2 Step 5: Final Form Rearranging: ∂ 2u 2 2∂ u = c ∂t2 ∂x2 q T where wave speed c = ρ Physical Interpretation Wave Speed q T The wave speed c = ρ depends on: Tension T : Higher tension → faster waves Linear density ρ: Higher density → slower waves Solutions General solution has the form: u(x, t) = f (x − ct) + g(x + ct) representing: f (x − ct): Waves moving right g(x + ct): Waves moving left 15 Key Properties Second order in both space and time Hyperbolic PDE Conserves energy Solutions maintain their shape while traveling Superposition principle applies (for linear wave equation) Applications The wave equation describes many physical phenomena: Vibrating strings (musical instruments) Sound waves in air Electromagnetic waves Water waves (in certain limits) Seismic waves Boundary Value Problems (BVPs) Definition A boundary value problem consists of: A differential equation (ordinary or partial) A domain where the solution is sought Additional conditions specified at the domain boundaries Types of Boundary Conditions Dirichlet Conditions The value of the solution is specified at the boundary: u(a) = α, u(b) = β Example: Temperature specified at the ends of a rod u(0) = T1 , u(L) = T2 Neumann Conditions The derivative of the solution is specified at the boundary: du du (a) = α, (b) = β dx dx Example: Insulated end of a rod (no heat flux) ∂u (0) = 0 ∂x 16 Mixed/Robin Conditions A combination of the solution and its derivative: du au(0) + b (0) = c dx Example: Convective cooling at boundary ∂u k (L) + h[u(L) − T∞ ] = 0 ∂x Initial Conditions Definition Initial conditions specify the state of a system at the starting time (usually t = 0) in time- dependent problems. Examples for Different PDEs Heat Equation For the heat equation: ∂u ∂ 2u =α 2 ∂t ∂x Initial condition (initial temperature distribution): u(x, 0) = f (x) Wave Equation For the wave equation: ∂ 2u 2 2∂ u = c ∂t2 ∂x2 Two initial conditions required: u(x, 0) = f (x) (initial position) ∂u (x, 0) = g(x) (initial velocity) ∂t Complete Initial-Boundary Value Problems Heat Conduction Example Consider heat conduction in a rod of length L: ∂u ∂ 2u = α 2 , 0 < x < L, t > 0 ∂t ∂x u(0, t) = T1 , u(L, t) = T2 (boundary conditions) u(x, 0) = f (x) (initial condition) 17 Wave Propagation Example Consider a vibrating string: ∂ 2u 2 2∂ u = c , 0 < x < L, t > 0 ∂t2 ∂x2 u(0, t) = 0, u(L, t) = 0 (fixed ends) u(x, 0) = f (x) (initial displacement) ∂u (x, 0) = g(x) (initial velocity) ∂t Importance in Applications Physical Significance Boundary conditions represent physical constraints or interactions at boundaries Initial conditions represent how system is prepared or starts Together they ensure unique solutions Connect mathematical models to measurable quantities Solution Requirements First-order equations in time need one initial condition Second-order equations in time need two initial conditions Number of boundary conditions depends on spatial order of equation Conditions must be compatible at corners of space-time domain Important Properties and Relationships Superposition Principle For linear, homogeneous PDEs: If u1 and u2 are solutions, then c1 u1 + c2 u2 is also a solution Allows decomposition of complex problems into simpler ones Boundary Conditions Number of required boundary conditions depends on order First-order: One boundary condition Second-order: Two boundary conditions nth-order: n boundary conditions 18 Solution Methods Linear PDEs: Wider range of analytical methods available – Separation of variables – Fourier transforms – Green’s functions Nonlinear PDEs: Often require numerical methods – Finite difference methods – Finite element methods – Perturbation methods Heat Equation - Separation of Variables Consider the heat equation on a finite rod [0, L]: ∂u ∂ 2u = α 2 , 0 < x < L, t > 0 ∂t ∂x u(0, t) = 0, u(L, t) = 0 (boundary conditions) u(x, 0) = f (x) (initial condition) Solution Steps Step 1: Assume Separated Solution Assume solution has the form: u(x, t) = X(x)T (t) Step 2: Substitute into PDE dT d2 X X(x) = αT (t) 2 dt dx Step 3: Separate Variables Divide both sides by αX(x)T (t): 1 dT 1 d2 X = = −λ αT dt X dx2 where −λ is the separation constant. Step 4: Solve Two ODEs Time equation: dT = −αλT =⇒ T (t) = ce−αλt dt Space equation: d2 X + λX = 0 dx2 19 with boundary conditions: X(0) = 0, X(L) = 0 Step 5: Solve Spatial Eigenvalue Problem The spatial equation with boundary conditions gives: nπx nπ 2 Xn (x) = sin , λn = L L where n = 1, 2, 3,... Step 6: Write General Solution The general solution is: ∞ nπx 2t X u(x, t) = cn sin e−α(nπ/L) n=1 L Finding Coefficients Step 7: Apply Initial Condition At t = 0: ∞ X nπx f (x) = u(x, 0) = cn sin n=1 L Step 8: Use Fourier Series The coefficients are found using orthogonality: 2 L Z nπx cn = f (x) sin dx L 0 L Final Solution The complete solution is: ∞ Z L X 2 nπx nπx 2t u(x, t) = f (x) sin dx sin e−α(nπ/L) n=1 L 0 L L Physical Interpretation Spatial Modes Each term in the series represents a mode: sin(nπx/L) gives spatial shape Higher modes (n large) oscillate more rapidly 20 Temporal Behavior Each mode decays exponentially with time Higher modes decay faster (due to n2 in exponent) System approaches steady state as t → ∞ Additional Resources Please take a look at these for a detailed workthrough: https://web.uvic.ca/~tbazett/diffyqs/heateq_section.html https://web.uvic.ca/~tbazett/diffyqs/we_section.html Solving PDEs using Fourier Transforms The Fourier Transform (FT) is a powerful technique for solving Partial Differential Equations (PDEs) by converting complex differential operations into simpler algebraic ones. Basic Principles The Fourier Transform converts derivatives into multiplication by powers of iω: ∂f F = iωF{f } ∂x 2 ∂ f F 2 = −ω 2 F{f } ∂x General Procedure 1. Take the Fourier Transform of both sides of the PDE 2. Solve the resulting algebraic equation 3. Take the inverse Fourier Transform to get the solution Example: Heat Equation Let’s solve the heat equation in one dimension: ∂u ∂ 2u =α 2 (5) ∂t ∂x Step 1: Apply Fourier Transform Let û(ω, t) = F{u(x, t)}. Taking the Fourier Transform with respect to x: ∂ û = −αω 2 û(ω, t) (6) ∂t 21 Step 2: Solve the ODE This is a first-order ODE with solution: 2t û(ω, t) = û(ω, 0)e−αω (7) where û(ω, 0) is the Fourier Transform of the initial condition. Step 3: Take Inverse Fourier Transform The solution in physical space is: 2 u(x, t) = F −1 {û(ω, 0)e−αω t } (8) For a specific initial condition u(x, 0) = f (x), the solution becomes: Z ∞ 1 (x−y)2 u(x, t) = √ f (y)e− 4αt dy (9) 4παt −∞ This is known as the fundamental solution or Green’s function solution of the heat equation. Solving PDEs using Laplace Transforms The Laplace Transform (LT) is particularly useful for solving PDEs with initial value problems. Unlike Fourier transforms which handle spatial derivatives, Laplace transforms are typically applied to the time variable. Basic Principles The Laplace Transform converts time derivatives into multiplication by powers of s: ∂f L = sL{f } − f (x, 0) ∂t 2 ∂ f ∂f L 2 = s2 L{f } − sf (x, 0) − (x, 0) ∂t ∂t General Procedure 1. Apply Laplace Transform with respect to time t 2. Solve the resulting ODE in terms of the spatial variable 3. Apply inverse Laplace Transform to obtain the solution Example: Wave Equation Consider the one-dimensional wave equation: ∂ 2u 2 2∂ u =c (10) ∂t2 ∂x2 with initial conditions: u(x, 0) = f (x) ∂u (x, 0) = g(x) ∂t 22 Step 1: Apply Laplace Transform Let U (x, s) = L{u(x, t)}. Taking the Laplace Transform: ∂ 2U s2 U (x, s) − sf (x) − g(x) = c2 (11) ∂x2 Step 2: Solve the ODE Rearranging: ∂ 2U s2 sf (x) + g(x) 2 − 2U = − (12) ∂x c c2 The general solution is: U (x, s) = A(s)esx/c + B(s)e−sx/c + particular solution (13) where A(s) and B(s) are determined from boundary conditions. Step 3: Take Inverse Laplace Transform For infinite domain with zero boundary conditions at infinity: 1 x+ct Z 1 u(x, t) = [f (x + ct) + f (x − ct)] + g(ξ)dξ (14) 2 2c x−ct This is D’Alembert’s solution to the wave equation. Advantages of Laplace Transform Particularly effective for problems with initial conditions Handles discontinuous and impulsive forcing functions well Transforms differential equations into algebraic equations Natural for problems with time-varying boundary conditions Common Applications 1. Heat conduction problems 2. Wave propagation 3. Beam deflection 4. Diffusion processes 5. Vibration analysis 23 Introduction Finite difference methods (FDM) are numerical techniques for solving partial differential equa- tions (PDEs) by approximating derivatives with differences between discrete points. The method involves: 1. Creating a grid or mesh over the domain 2. Replacing continuous derivatives with discrete approximations 3. Transforming the PDE into a system of algebraic equations Basic Concepts Consider the one-dimensional heat equation: ∂u ∂ 2u =α 2 (15) ∂t ∂x Spatial Discretization For the spatial second derivative, we use the central difference approximation: ∂ 2u u(x + h) − 2u(x) + u(x − h) 2 ≈ (16) ∂x h2 where h is the grid spacing. Temporal Discretization For the time derivative, we use a forward difference: ∂u u(t + ∆t) − u(t) ≈ (17) ∂t ∆t where ∆t is the time step. Numerical Implementation Combining these approximations, we get the explicit scheme: un+1 i = uni + r(uni+1 − 2uni + uni−1 ) (18) where r = α∆t/h2 is the stability parameter. Python Implementation import numpy as np def s o l v e h e a t e q u a t i o n ( i n i t i a l t e m p , dx , dt , alpha , t o t a l t i m e ) : # C a l c u l a t e number o f time s t e p s nt = int ( t o t a l t i m e / dt ) nx = len ( i n i t i a l t e m p ) 24 # I n i t i a l i z e s o l u t i o n array u = np. z e r o s ( ( nt , nx ) ) u = initial temp # S t a b i l i t y condition r = a lpha ∗ dt / ( dx ∗∗ 2) if r > 0.5: r a i s e Val ueErro r ( ” S o l u t i o n may be u n s t a b l e ” ) # S o l v e u s i n g e x p l i c i t scheme for t in range ( 0 , nt −1): for i in range ( 1 , nx −1): u [ t +1, i ] = u [ t , i ] + r ∗ ( u [ t , i +1] − 2∗u [ t , i ] + u [ t , i −1]) # Apply boundary c o n d i t i o n s u [ t +1 ,0] = u [ t +1 ,1] # Neumann a t x=0 u [ t +1 , −1] = u [ t +1 , −2] # Neumann a t x=L return u Important Considerations Accuracy The truncation error depends on the grid spacing: Spatial error: O(h2 ) for central difference Temporal error: O(∆t) for forward difference Stability The explicit scheme is conditionally stable. For the heat equation: α∆t 1 r= 2 ≤ (19) h 2 Boundary Conditions Common types include: Dirichlet: u(xb , t) = g(t) Neumann: ∂u ∂x (xb , t) = g(t) Periodic: u(x0 , t) = u(xN , t) 25 Alternative Schemes Implicit (Backward Difference) More stable but requires solving a system: −run+1 n+1 i−1 + (1 + 2r)ui − run+1 n i+1 = ui (20) Crank-Nicolson Semi-implicit scheme with better accuracy: n n+1 un+1 − uni ∂ 2u ∂ 2u i α α = + (21) ∆t 2 ∂x2 2 ∂x2 Explicit Scheme (FTCS) vs. Implicit Scheme (BTCS) In finite difference methods for solving partial differential equations (PDEs), two primary ap- proaches exist: explicit and implicit schemes. These differ fundamentally in how they compute the solution at the next time step. Explicit Scheme (FTCS) Basic Concept In explicit schemes, we compute the solution at time step n + 1 directly using known values from time step n. Consider the heat equation: ∂u ∂ 2u =α 2 (22) ∂t ∂x The explicit (forward-time) scheme takes the form: un+1 i = uni + r(uni+1 − 2uni + uni−1 ) (23) where r = α∆t/(∆x)2 is the stability parameter. Stencil Representation The explicit scheme uses this computational stencil: Advantages Simple to implement Computationally inexpensive per time step Parallelizes easily 26 FTCS Stencil Image from : https://web.cecs.pdx.edu/~gerry/class/ME448/notes/1Dmodels/finiteDifferenceHeatEquation/ Disadvantages Conditionally stable (r ≤ 1 2 for heat equation) Requires small time steps First-order accurate in time Implicit Schemes (BTCS) Basic Concept Implicit schemes compute the solution at time step n + 1 by solving a system of equations involving unknown values at the new time step. The backward-time scheme is: un+1 i − r(un+1 n+1 i+1 − 2ui + un+1 n i−1 ) = ui (24) This can be rewritten as: −run+1 n+1 i−1 + (1 + 2r)ui − run+1 n i+1 = ui (25) Matrix Form The system becomes:    n+1    1 + 2r −r 0 ··· 0 u1 un1  −r 1 + 2r −r ··· 0   n+1   n     u2   u2  ............. .  .   ..  = ..  (26)  .. −r  un+1     n   0 ··· −r 1 + 2r N −1  uN −1  n+1 0 ··· 0 −r 1 + 2r uN unN Advantages Unconditionally stable Allows larger time steps Better handling of stiff problems 27 BTCS Stencil Image from : https://web.cecs.pdx.edu/~gerry/class/ME448/notes/1Dmodels/finiteDifferenceHeatEquation/ Disadvantages More computationally expensive per time step Requires solving a system of equations More complex to implement Harder to parallelize Comparison of Stability The stability of these schemes can be analyzed using von Neumann stability analysis. For the heat equation: Explicit Scheme: |G(ξ)| = |1 − 4r sin2 (ξ∆x/2)| ≤ 1 (27) This leads to the stability condition: 1 r≤ (28) 2 Implicit Scheme: 1 |G(ξ)| = ≤1 (29) |1 + 4r sin2 (ξ∆x/2)| This is always satisfied for any positive r, making the scheme unconditionally stable. Practical Considerations Choice of Scheme The choice between explicit and implicit schemes depends on: Required accuracy Stability constraints Available computational resources Problem characteristics (stiffness) Time step requirements 28 Hybrid Approaches Crank-Nicolson scheme combines explicit and implicit approaches: n n+1 un+1 − uni α ∂ 2u α ∂ 2u i = + (30) ∆t 2 ∂x2 2 ∂x2 This provides: Second-order accuracy in time Unconditional stability Better balance of accuracy and stability Finite Element Methods (FEM) Finite Element Methods (FEM) are numerical techniques for approximating solutions to partial differential equations (PDEs) by dividing a complex domain into simpler subdomains called elements. These elements connect at points called nodes. Steps in FEM 1. Discretization: Breaking the complex geometry into smaller elements, creating a so called “mesh”. 2. Element Analysis: For each small element, we write equations that approximately describe the behavior we’re interested in (like heat transfer, stress, or fluid flow) 3. Assembly: Combining all these small element equations into one large system of equa- tions that represents the entire problem 4. Solving: Using numerical methods to solve this large system of equations 5. Post-processing: Interpreting the results and visualizing them A major advantage of FEM is its ability to handle irregular geometries and complex bound- ary conditions that would be extremely difficult or impossible to solve analytically. It’s widely used in structural analysis, heat transfer, fluid dynamics, and electromagnetic problems. Weak Form of the PDE The weak form of a Partial Differential Equation (PDE) is a reformulation of the original “strong form” that relaxes the continuity requirements on the solution and is particularly suitable for finite element analysis. Strong Form vs Weak Form The key differences between the strong and weak forms are: The strong form requires the solution to satisfy the PDE at every point in the domain The weak form requires the solution to satisfy the equation in an average sense when multiplied by test functions and integrated over the domain 29 Derivation Process The weak form is derived through the following steps: 1. Multiply the PDE by a test function (typically denoted as v) 2. Integrate over the domain 3. Use integration by parts to reduce the order of derivatives 4. Apply boundary conditions Example: Poisson Equation Consider a simple 1D Poisson equation: Strong Form d2 u − = f in Ω dx2 u = 0 on ∂Ω Deriving the Weak Form 1. Multiply by test function v and integrate: d2 u Z Z − 2 v dx = f v dx Ω dx Ω 2. Apply integration by parts: Z 1 Z du dv du dx − v = f v dx Ω dx dx dx 0 Ω 3. With zero boundary conditions, the boundary term vanishes, giving the weak form: Z Z du dv dx = f v dx Ω dx dx Ω Advantages of Weak Form The weak form offers several advantages: Lower continuity requirements on the solution Natural incorporation of boundary conditions Mathematical framework better suited for existence and uniqueness proofs Direct connection to variational principles Forms the basis for Galerkin finite element methods 30 Mathematical Framework The weak form transforms the problem into finding u ∈ H 1 (Ω) (a Sobolev space) such that: a(u, v) = l(v) ∀v ∈ H01 (Ω) where: Z du dv a(u, v) = dx dx dx ZΩ l(v) = f v dx Ω This formulation requires less smoothness than the strong form, which needs u ∈ C 2 (Ω). Additional Material https://www.comsol.com/blogs/brief-introduction-weak-form https://www.simscale.com/forum/t/the-finite-element-method-fundamentals-strong-and-weak-form-for-1d-problems-5/57940 Variational Formulation of the Poisson Equation 1. The Poisson Equation The Poisson equation in a domain Ω ⊂ Rd is given by: −∇ · (∇u) = f in Ω, with boundary conditions: Dirichlet boundary condition: u = uD on ∂ΩD , Neumann boundary condition: ∂u ∂n = g on ∂ΩN , where ∂Ω = ∂ΩD ∪ ∂ΩN , and n is the outward unit normal vector. Here: u(x): Unknown solution, f (x): Source term, ∇u = ∂x ∂u 1 ,... , ∂u ∂xd : Gradient. 2. Variational Formulation To derive the variational formulation: 1. Multiply by a test function v: Multiply both sides of the equation by a test function v ∈ V , where V is a space of sufficiently smooth functions (e.g., H 1 (Ω)): Z Z − (∇ · ∇u)v dΩ = f v dΩ. Ω Ω 31 2. Apply integration by parts: Using the divergence theorem, transform the second- order derivative term into a first-order term: Z Z Z ∇u · ∇v dΩ − v(∇u · n) dS = f v dΩ. Ω ∂Ω Ω Here: ∂Ω : Integral over the boundary of Ω, R ∇u · n: Flux term. 3. Incorporate boundary conditions: For Neumann boundary conditions (∇u · n = g on ∂ΩN ), substitute into the boundary integral: Z Z Z ∇u · ∇v dΩ = f v dΩ + gv dS. Ω Ω ∂ΩN For Dirichlet boundary conditions (u = uD on ∂ΩD ), restrict the solution space V to functions that vanish on ∂ΩD : V = {v ∈ H 1 (Ω) | v|∂ΩD = 0}. 3. Weak Formulation The weak (variational) formulation is: Find u ∈ V such that: Z Z Z ∇u · ∇v dΩ = f v dΩ + gv dS ∀v ∈ V. Ω Ω ∂ΩN Here: u ∈ V : The solution u satisfies the Dirichlet boundary conditions. v ∈ V : The test functions v vanish on ∂ΩD. 4. Finite Element Discretization To solve numerically: 1. Discretize the domain: Divide Ω into smaller elements (e.g., triangles or tetrahedra). 2. Choose a finite-dimensional subspace: Select Vh ⊂ V , a space of piecewise polyno- mial functions (e.g., linear or quadratic basis functions). 3. Approximate the solution and test functions: Let u ≈ uh ∈ Vh and v ≈ vh ∈ Vh. 4. Solve the resulting linear system: Au = b, where: A: Stiffness matrix from Ω ∇uh · ∇vh dΩ, R b: Load vector from Ω f vh dΩ + ∂ΩN gvh dS, R R u: Vector of nodal values of the approximate solution. 32 Weighted Residual Method (WRM) The Weighted Residual Method (WRM) is a fundamental mathematical technique used in the Finite Element Method (FEM) for finding approximate solutions to differential equations. This document explains the core concepts and implementation of WRM in FEM analysis. Basic Concept Consider a differential equation in operator form: L(u) = f (31) where L is a differential operator. When we seek an approximate solution uh , we obtain a residual: R = L(uh ) − f ̸= 0 (32) In the weighted residual method, we begin with a differential equation and its boundary conditions. The method involves: Assuming an approximate solution with unknown parameters Recognizing that this approximation creates a residual Minimizing this residual using weight functions Mathematical Formulation Approximation We approximate the solution using trial functions (basis functions): n X u(x) ≈ uh (x) = Ni (x)ai (33) i=1 where: Ni (x) are the basis functions ai are unknown coefficients n is the number of basis functions Residual The residual R(x) is defined as: R(x) = L(uh ) − f (34) where: L is the differential operator f is the known function 33 Weighted Integral Form The weighted residual equation is: Z wi (x)R(x)dx = 0 (35) Ω where wi (x) are the weight functions and Ω is the domain. Common Weight Functions Galerkin Method The most commonly used method in FEM: wi (x) = Ni (x) (36) Collocation Method Uses Dirac delta functions: wi (x) = δ(x − xi ) (37) Least Squares Method Uses derivatives of the residual: ∂R wi (x) = (38) ∂ai Subdomain Method Uses step functions over subdomains: ( 1 in subdomain i wi (x) = (39) 0 elsewhere Advantages Provides systematic approach to FEM formulation Handles complex geometries and boundary conditions Offers mathematical rigor and consistency Often leads to symmetric matrices (particularly with Galerkin’s method) Application in FEM The weighted residual method in FEM: Transforms continuous problems into discrete systems Creates solvable systems of algebraic equations Provides framework for error analysis Enables convergence studies 34 Conclusion The Galerkin method, as a special case of WRM, is particularly valuable in FEM due to its optimal convergence properties and the symmetric matrices it produces. The success of the method depends significantly on the choice of appropriate basis and weight functions that satisfy the required mathematical conditions. The Galerkin Method Fundamental Principle The Galerkin method enforces orthogonality between the residual and all test functions in the approximation space: Z R(x)vi (x) dx = 0 for all test functions vi (40) Ω Implementation 1. Approximate the solution using basis functions: n X uh (x) = cj ϕj (x) (41) j=1 2. Form the residual equation: n ! X R(x) = L cj ϕj (x) − f (x) (42) j=1 3. Apply the Galerkin condition: Z " Xn ! # L cj ϕj (x) − f (x) ϕi (x) dx = 0 for i = 1, 2,... , n (43) Ω j=1 Connection to Variational Form Theorem 1 For self-adjoint operators, the Galerkin method produces equations identical to those derived from variational principles. Example: Poisson Equation Consider −u′′ (x) = f (x): 1. Residual: ! n X R=− cj ϕ′′j (x) − f (x) (44) j=1 2. Galerkin condition: Z " n X ! # − cj ϕ′′j (x) − f (x) ϕi (x) dx = 0 (45) Ω j=1 35 3. After integration by parts: Z n ! Z X cj ϕ′j (x) ϕ′i (x) dx = f (x)ϕi (x) dx (46) Ω j=1 Ω Properties and Advantages Generality: Applicable to non-variational problems Optimality: Produces optimal approximations in energy norm for symmetric problems Matrix Properties: Results in symmetric matrices for self-adjoint operators Boundary Conditions: Natural treatment of essential and natural boundary conditions Relationship to Other Methods Definition 2 (General Weighted Residual Method) A weighted residual method seeks to minimize the residual in some sense by choosing appropriate weight functions wi (x): Z R(x)wi (x) dx = 0 (47) Ω Remark 3 The Galerkin method is distinguished by its choice of weight functions: wi (x) = ϕi (x) (48) where ϕi (x) are the same functions used in the trial space. Error Analysis For symmetric, positive-definite problems, the Galerkin method produces approximations that are optimal in the energy norm: ∥u − uh ∥E = min ∥u − vh ∥E (49) vh ∈Vh where ∥ · ∥E denotes the energy norm associated with the bilinear form. 36

Introduction to Mathematical Modeling and Partial Differential Equations PDF

Document Details

Tags

Related

Summary

Full Transcript

Upgrade to continue