Theoretical Physics Project 4- Partial differential equations (PDEs) Solution

Project 4: Partial differential equations (PDEs)
1. Laplace equation (10 points) 100 V
Consider the two-dimensional setup of a square box with a side length of L, where the top edge is set to a potential of 100V and the other three edges are grounded, i.e. at a potential of zero. Your task is to find the potential within the box, where it solves the Laplace equation
∆ϕ(x,y) = 0. (1) a) Write a program that computes the potential ϕ(x,y)
within the box with a discretisation length of ∆x = L/100. Implement the Jacobi, Gauß–Seidel and SOR methods. As a stopping criterion, assert that the error of the discretised Laplace equation is smaller than ϵmax = 10−3 V everywhere. The error is defined as ϵi,j = |ϕi,j −(ϕi+1,j +ϕi−1,j +ϕi,j+1 +ϕi,j−1)/4|. Plot the average and maximum error versus the number of iterations for all algorithms. For SOR, use four different over-relaxation parameter values α = 0.5,1.0,1.25,1.5,1.75 and 1.99.
Check if it still converges for α ≥ 2.0. Discuss your results. (6 points)
b) The solution of this problem can be written as an infinite series:
. (2)
Plot this solution for 1, 10, 100 and 1000 terms. For terms with large n, where the sinh terms become large and might cause overflow, you might want to use the fact that the ratio of sinh terms approaches exp(−nπx/L) with n → ∞. Plot the difference between one of the iterative solutions and the “infinite” series solution using 1000 terms. Discuss your results. (4 points)
2. Diffusion (16 points)
Consider a metallic rod of a finite length L and a small radius, which is isolated at its side, but not at its ends, where it is placed in contact with ice water at 0◦C:

To simulate the temperature flow, we will need to solve the diffusion PDE
, (3)
with the thermal conductivity K, the heat capacity C and the density ρ. Use L = 1, K = 210, C = 900 and ρ = 2700, and, if not otherwise specified, ∆x = 0.01 and ∆t = 0.1 with a number of time steps Nt = 10000. The rod’s temperature distribution is initially set to T(x,0) = sin(πx/L).
a) Simulate the system using the FTCS algorithm and plot T(x,t). (8 points)
b) The analytical solution is given by
πx Å π2Kt
. (4)
With that, we can define our simulation error as
Nx−1
, (5)
j=1
where Nx = L/∆x. Study the behaviour of ϵ(t = 100) varying ∆t between 0.001 and 0.7 while keeping ∆x fixed. What happens? Explain your finding. (4 points)
c) Now, implement the implicit Euler backward, the Crank–Nicolson and Dufort– Frankel algorithms and repeat the analysis for ϵ(t = 100) for all of them. Discuss your results comparing the four algorithms, in particular with respect to the scaling behaviour of the error with ∆t. (4 points)
3. Solitons (16 points)
Water waves in shallow, narrow channels can be described by the Koorteweg–de Vries (KdeV) equation
. (6)
The non-linear term leads to a sharpening of the wave and ultimately to a shock wave. In contrast, the ∂3u/∂x3 term produces broadening. For the proper parameters and initial conditions, the two effects exactly balance each other, and a stable wave is formed, which is called a “soliton”. These stable solitons almost behave as particles, and appear in many areas of physics. For more details on the numerical discovery of this fascinating phenomenon, see “Computer Simulations Led to Discovery of Solitons”. For a deeper dive, you can also download the original paper from Stud.IP ▶ Ubung: Methods o¨ f Computational Physics ▶ Files ▶ Additional Material.
Inserting the traveling wave ansatz u(x,t) = u(x − ct) gives a solvable ODE with the (non-trivial) solution
, (7)
where ξ = x−ct is the phase. Note that the amplitude is proportional to the propagation speed c. Discretising the KdeV equation gives the algorithm
(8)
.
To calculate u2j, the first term on the right-hand side is not yet available, so we use a simpler forward difference scheme for the time step, which gives
(9)
.
In addition to the first time step, we also need to treat and unNmax separately, since also they can not be updated using Eq. (8). Given the initial state we will use below, a simple method is to use the constant boundary values = 1 and unNmax = 0. This leaves and unNmax−1. These are still to be calculated using Eqs. (8) and (9), but with the replacements (when calculating ) and unNmax+1 → unNmax (when calculating unNmax−1) on the right-hand sides.
a) Derive Eq. (8). For the first derivates, use the central difference approximations
∂u/∂t ≈ (unj +1 − unj −1)/(2∆t) and ). To find an approximation for the third derivative expand u(x,t) to including O(∆x4) about the four points u(x ± ∆x,t) and u(x ± 2∆x,t), then solve for ∂3u/∂x3. Finally, for the non-differentiated u(x,t) in the second term of Eq. (8), use the average over three adjacent points, i.e. 3. In your derivation, keep track of the truncation error, to prove that Eq. (8) is of order O(∆t3) + O(∆t∆x2). (4 points)
b) Implement the above algorithm solving the KdeV equation for the initial condition
1 ï Åx − 25ãò
x,t . (10)
The parameters are given by ϵ = 0.2 and µ = 0.1 and the step sizes by ∆x = 0.4 and ∆t = 0.1. Use 130 steps in x, such that our simulation volume is given by L = 130∆x = 52. Verify that these constants satisfy the stability condition
1. (11)
∆x (∆x)2
Store the solution every 250 or so time steps and run the simulation for about 2000 time steps. Plot the disturbance u versus position and versus time. Scott Russell observed in 1834 in the Edinburgh-Glasgow canal that an initial, arbitrary wave form set in motion evolves into two or more solitary waves that move at different velocities. Can you confirm his observation? Into how many solitary waves does your initial wave form break up into? (8 points)
c) Explore what happens when a tall soliton collides with a short one. Do they bounce off each other? Do they go through each other? Do they interfere? Do they destroy each other? Does the tall soliton still move faster than the short one after collision? Start off by placing a tall soliton of height 0.8 at x = 12, and a small soliton in front of it at x = 26:
. (12)
Make sure to adapt the constant boundary values for and unNmax for this simulation. (4 points)