Waves On A String Simulation
A very wide range of physical processes pb to wave movement, where signals are propagated through a medium in infinite and time, normally with little or no permanent movement of the medium itself. The shape of the signals may undergo changes as they travel through matter, but ordinarily not so much that the signals cannot exist recognized at some later point in space and time. Many types of wave move can be described by the equation \(u_{tt}=\nabla\cdot (c^2\nabla u) + f\), which we will solve in the forthcoming text by finite difference methods.
Simulation of waves on a string¶
We begin our report of wave equations past simulating one-dimensional waves on a cord, say on a guitar or violin. Allow the string in the undeformed state coincide with the interval \([0,L]\) on the \(ten\) axis, and let \(u(ten,t)\) be the displacement at time \(t\) in the \(y\) management of a point initially at \(x\). The displacement function \(u\) is governed by the mathematical model
\[ \begin{equation} \frac{\partial^two u}{\partial t^2} = c^ii \frac{\partial^2 u}{\partial x^2}, \quad x\in (0,Fifty),\ t\in (0,T] \characterization{wave:pde1} \tag{1} \end{equation} \]
\[ \brainstorm{equation} u(x,0) = I(x), \quad 10\in [0,50] \label{wave:pde1:ic:u} \tag{2} \end{equation} \]
\[ \begin{equation} \frac{\partial}{\partial t}u(x,0) = 0, \quad x\in [0,50] \label{moving ridge:pde1:ic:ut} \tag{iii} \end{equation} \]
\[ \begin{equation} u(0,t) = 0, \quad t\in (0,T] \characterization{wave:pde1:bc:0} \tag{4} \finish{equation} \]
\[ \begin{equation} u(50,t) = 0, \quad t\in (0,T] \label{moving ridge:pde1:bc:L} \tag{5} \end{equation} \]
The abiding \(c\) and the function \(I(x)\) must exist prescribed.
Equation (i) is known as the 1-dimensional wave equation. Since this PDE contains a second-order derivative in time, nosotros need ii initial weather condition. The condition (ii) specifies the initial shape of the string, \(I(x)\), and (three) expresses that the initial velocity of the string is naught. In addition, PDEs need purlieus conditions, given here as (iv) and (5). These two atmospheric condition specify that the string is fixed at the ends, i.e., that the deportation \(u\) is nil.
The solution \(u(ten,t)\) varies in infinite and time and describes waves that motion with velocity \(c\) to the left and right.
Sometimes we will use a more than compact note for the partial derivatives to save space:
\[ \begin{equation} u_t = \frac{\fractional u}{\fractional t}, \quad u_{tt} = \frac{\partial^ii u}{\partial t^ii}, \label{_auto1} \tag{six} \end{equation} \]
and similar expressions for derivatives with respect to other variables. So the wave equation tin be written compactly as \(u_{tt} = c^2u_{xx}\).
The PDE problem (1)-(five) will at present be discretized in space and time by a finite difference method.
Discretizing the domain¶
The temporal domain \([0,T]\) is represented by a finite number of mesh points
\[ \begin{equation} 0 = t_0 < t_1 < t_2 < \cdots < t_{N_t-1} < t_{N_t} = T \characterization{_auto2} \tag{7} \end{equation} \]
Similarly, the spatial domain \([0,L]\) is replaced by a set of mesh points
\[ \begin{equation} 0 = x_0 < x_1 < x_2 < \cdots < x_{N_x-1} < x_{N_x} = L \label{_auto3} \tag{viii} \terminate{equation} \]
Ane may view the mesh equally two-dimensional in the \(x,t\) plane, consisting of points \((x_i, t_n)\), with \(i=0,\ldots,N_x\) and \(n=0,\ldots,N_t\).
Uniform meshes¶
For uniformly distributed mesh points we can introduce the constant mesh spacings \(\Delta t\) and \(\Delta x\). We accept that
\[ \begin{equation} x_i = i\Delta x,\ i=0,\ldots,N_x,\quad t_n = n\Delta t,\ n=0,\ldots,N_t \label{_auto4} \tag{9} \end{equation} \]
We also have that \(\Delta x = x_i-x_{i-i}\), \(i=1,\ldots,N_x\), and \(\Delta t = t_n - t_{n-1}\), \(n=one,\ldots,N_t\). Figure displays a mesh in the \(x,t\) plane with \(N_t=five\), \(N_x=5\), and constant mesh spacings.
The discrete solution¶
The solution \(u(x,t)\) is sought at the mesh points. We introduce the mesh function \(u_i^n\), which approximates the verbal solution at the mesh point \((x_i,t_n)\) for \(i=0,\ldots,N_x\) and \(n=0,\ldots,N_t\). Using the finite divergence method, we shall develop algebraic equations for computing the mesh role.
Fulfilling the equation at the mesh points¶
In the finite difference method, nosotros relax the status that (1) holds at all points in the space-time domain \((0,50)\times (0,T]\) to the requirement that the PDE is fulfilled at the interior mesh points but:
\[ \brainstorm{equation} \frac{\partial^2}{\partial t^2} u(x_i, t_n) = c^ii\frac{\partial^2}{\partial x^two} u(x_i, t_n), \label{wave:pde1:step2} \tag{10} \end{equation} \]
for \(i=1,\ldots,N_x-ane\) and \(due north=1,\ldots,N_t-1\). For \(n=0\) we have the initial conditions \(u=I(x)\) and \(u_t=0\), and at the boundaries \(i=0,N_x\) we have the boundary condition \(u=0\).
Replacing derivatives by finite differences¶
The second-order derivatives can be replaced by key differences. The most widely used difference approximation of the 2nd-order derivative is
\[ \frac{\fractional^ii}{\partial t^2}u(x_i,t_n)\approx \frac{u_i^{n+one} - 2u_i^northward + u^{n-ane}_i}{\Delta t^2}\ \]
mathcal{I}_t is convenient to introduce the finite difference operator annotation
\[ [D_tD_t u]^n_i = \frac{u_i^{due north+1} - 2u_i^n + u^{n-1}_i}{\Delta t^2} \]
A like approximation of the 2d-order derivative in the \(10\) direction reads
\[ \frac{\partial^2}{\partial x^2}u(x_i,t_n)\approx \frac{u_{i+ane}^{n} - 2u_i^n + u^{n}_{i-1}}{\Delta x^ii} = [D_xD_x u]^n_i \]
Algebraic version of the PDE¶
We can now replace the derivatives in (10) and get
\[ \begin{equation} \frac{u_i^{n+1} - 2u_i^n + u^{n-1}_i}{\Delta t^2} = c^2\frac{u_{i+1}^{north} - 2u_i^n + u^{n}_{i-i}}{\Delta ten^two}, \label{wave:pde1:step3b} \tag{11} \end{equation} \]
or written more than compactly using the operator notation:
\[ \begin{equation} [D_tD_t u = c^2 D_xD_x]^{due north}_i \label{wave:pde1:step3a} \tag{12} \end{equation} \]
Interpretation of the equation as a stencil¶
A characteristic feature of (11) is that it involves \(u\) values from neighboring points just: \(u_i^{n+1}\), \(u^n_{i\pm ane}\), \(u^n_i\), and \(u^{n-1}_i\). The circles in Effigy illustrate such neighboring mesh points that contribute to an algebraic equation. In this particular case, we have sampled the PDE at the point \((2,two)\) and constructed (11), which and so involves a coupling of \(u_1^2\), \(u_2^3\), \(u_2^ii\), \(u_2^1\), and \(u_3^2\). The term stencil is often used about the algebraic equation at a mesh point, and the geometry of a typical stencil is illustrated in Figure. Ane likewise often refers to the algebraic equations every bit detached equations, (finite) difference equations or a finite difference scheme.
Mesh in space and time. The circles bear witness points connected in a finite difference equation.
Algebraic version of the initial conditions¶
We also need to replace the derivative in the initial condition (3) by a finite deviation approximation. A centered deviation of the type
\[ \frac{\fractional}{\partial t} u(x_i,t_0)\approx \frac{u^1_i - u^{-1}_i}{2\Delta t} = [D_{2t} u]^0_i, \]
seems appropriate. Writing out this equation and ordering the terms requite
\[ \begin{equation} u^{-1}_i=u^{1}_i,\quad i=0,\ldots,N_x \label{wave:pde1:step3c} \tag{13} \end{equation} \]
The other initial condition can be computed by
\[ u_i^0 = I(x_i),\quad i=0,\ldots,N_x \]
Formulating a recursive algorithm¶
We assume that \(u^n_i\) and \(u^{north-ane}_i\) are available for \(i=0,\ldots,N_x\). The merely unknown quantity in (xi) is therefore \(u^{n+one}_i\), which we now can solve for:
\[ \begin{equation} u^{n+i}_i = -u^{due north-1}_i + 2u^n_i + C^2 \left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-one}\correct) \label{wave:pde1:step4} \tag{14} \end{equation} \]
Nosotros have hither introduced the parameter
\[ \begin{equation} C = c\frac{\Delta t}{\Delta x}, \label{_auto5} \tag{15} \finish{equation} \]
known every bit the Courant number.
\(C\) is the fundamental parameter in the discrete moving ridge equation.
We encounter that the discrete version of the PDE features simply one parameter, \(C\), which is therefore the key parameter, together with \(N_x\), that governs the quality of the numerical solution (come across the section Analysis of the difference equations for details). Both the primary physical parameter \(c\) and the numerical parameters \(\Delta ten\) and \(\Delta t\) are lumped together in \(C\). Note that \(C\) is a dimensionless parameter.
Given that \(u^{n-1}_i\) and \(u^n_i\) are known for \(i=0,\ldots,N_x\), nosotros find new values at the next time level by applying the formula (14) for \(i=1,\ldots,N_x-i\). Effigy illustrates the points that are used to compute \(u^3_2\). For the boundary points, \(i=0\) and \(i=N_x\), we apply the boundary conditions \(u_i^{n+one}=0\).
Even though sound reasoning leads upward to (14), in that location is still a pocket-size challenge with information technology that needs to be resolved. Recall of the very beginning computational stride to exist made. The scheme (fourteen) is supposed to starting time at \(n=1\), which ways that we compute \(u^two\) from \(u^i\) and \(u^0\). Unfortunately, we practise not know the value of \(u^1\), so how to proceed? A standard procedure in such cases is to apply (14) also for \(north=0\). This immediately seems strange, since it involves \(u^{-1}_i\), which is an undefined quantity outside the time mesh (and the time domain). However, we can use the initial condition (thirteen) in combination with (14) when \(northward=0\) to eliminate \(u^{-one}_i\) and arrive at a special formula for \(u_i^1\):
\[ \begin{equation} u_i^1 = u^0_i - \frac{1}{2} C^two\left(u^{0}_{i+1}-2u^{0}_{i} + u^{0}_{i-i}\right) \label{moving ridge:pde1:step4:1} \tag{16} \cease{equation} \]
Figure illustrates how (xvi) connects four instead of 5 points: \(u^1_2\), \(u_1^0\), \(u_2^0\), and \(u_3^0\).
Modified stencil for the get-go time step.
We can now summarize the computational algorithm:
-
Compute \(u^0_i=I(x_i)\) for \(i=0,\ldots,N_x\)
-
Compute \(u^1_i\) by (16) for \(i=i,2,\ldots,N_x-one\) and prepare \(u_i^ane=0\) for the boundary points given by \(i=0\) and \(i=N_x\),
-
For each time level \(due north=1,two,\ldots,N_t-1\)
a. apply (14) to notice \(u^{n+1}_i\) for \(i=one,\ldots,N_x-one\)
b. set \(u^{north+1}_i=0\) for the boundary points having \(i=0\), \(i=N_x\).
The algorithm essentially consists of moving a finite divergence stencil through all the mesh points, which tin can be seen every bit an blitheness in a web page or a moving picture file.
Sketch of an implementation¶
We start past defining some constants that will be used throughout our Devito lawmaking.
import numpy as np # Given mesh points as arrays ten and t (10[i], t[n]), # constant c and function I for initial status x = np . linspace ( 0 , 2 , 101 ) t = np . linspace ( 0 , 2 , 101 ) c = 1 I = lambda 10 : np . sin ( x ) dx = x [ 1 ] - x [ 0 ] dt = t [ ane ] - t [ 0 ] C = c * dt / dx # Courant number Nx = len ( x ) - 1 Nt = len ( t ) - 1 C2 = C ** 2 # Help variable in the scheme Fifty = 2.
Next, nosotros define our 1D computational filigree and create a function u every bit a symbolic devito.TimeFunction . We need to specify the space_order as ii since our wave equation involves second-order derivatives with respect to \(10\). Similarly, nosotros specify the time_order as two, equally our equation involves second-order derivatives with respect to \(t\). Setting these parameters allows usa to use u.dx2 and u.dt2 .
from devito import Grid , TimeFunction # Initialise `u` for space and time order 2, using initialisation function I grid = Filigree ( shape = ( Nx + 1 ), extent = ( L )) u = TimeFunction ( proper noun = 'u' , grid = grid , time_order = 2 , space_order = ii ) u . data [:,:] = I ( x [:])
Now that we have initialised u , we can solve our moving ridge equation for the unknown quantity \(u^{n+1}_i\) using forward and backward differences in infinite and time.
from devito import Constant , Eq , solve # Prepare wave equation and solve for forrad stencil bespeak in fourth dimension pde = ( 1 / c ** ii ) * u . dt2 - u . dx2 stencil = Eq ( u . forward , solve ( pde , u . forward )) print ( "LHS: %south " % stencil . lhs ) impress ( "RHS: %s " % stencil . rhs )
LHS: u(t + dt, x) RHS: 1.0*dt**two*(-2.0*u(t, x)/h_x**ii + u(t, x - h_x)/h_x**two + u(t, x + h_x)/h_x**two + 2.0*u(t, x)/dt**2 - 1.0*u(t - dt, 10)/dt**2)
Great! From these print statements, nosotros tin encounter that Devito has taken the moving ridge equation in (1) and solved information technology for \(u^{due north+i}_i\), giving us equation (xiv). Notation that dx is denoted equally h_x , while u(t, x) , u(t, x - h_x) and u(t, x + h_x) denote the equivalent of \(u^{due north}_{i}\), \(u^{due north}_{i-1}\) and \(u^{due north}_{i+1}\) respectively.
We as well need to create a separate stencil for the outset timestep, where nosotros substitute \(u^{1}_i\) for \(u^{-ane}_i\), as given in (13).
stencil_init = stencil . subs ( u . astern , u . forrard )
Now nosotros can create expressions for our purlieus weather condition and build the operator. The results are plotted below.
#NBVAL_IGNORE_OUTPUT from devito import Operator t_s = grid . stepping_dim # Purlieus conditions bc = [ Eq ( u [ t_s + one , 0 ], 0 )] bc += [ Eq ( u [ t_s + 1 , Nx ], 0 )] # Defining one Operator for initial timestep and i for the rest op_init = Operator ([ stencil_init ] + bc ) op = Operator ([ stencil ] + bc ) op_init . apply ( time_M = 1 , dt = dt ) op . employ ( time_m = one , time_M = Nt , dt = dt )
Data blazon float64 of runtime value `dt` does non lucifer the Constant information type <class 'numpy.float32'> Operator `Kernel` run in 0.01 s Information blazon float64 of runtime value `dt` does non match the Constant data type <form 'numpy.float32'> Operator `Kernel` run in 0.01 southward PerformanceSummary([(PerfKey(name='section0', rank=None), PerfEntry(time=two.8000000000000003e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))]) We tin plot our results using matplotlib :
import matplotlib.pyplot every bit plt plt . plot ( x , u . data [ - i ]) plt . xlabel ( 'x' ) plt . ylabel ( 'u' ) plt . show ()
Verification¶
Before implementing the algorithm, it is convenient to add together a source term to the PDE (1), since that gives united states of america more than freedom in finding test problems for verification. Physically, a source term acts as a generator for waves in the interior of the domain.
A slightly generalized model problem¶
Nosotros now address the following extended initial-boundary value problem for one-dimensional wave phenomena:
\[ \brainstorm{equation} u_{tt} = c^2 u_{twenty} + f(x,t), \quad x\in (0,L),\ t\in (0,T] \characterization{wave:pde2} \tag{17} \cease{equation} \]
\[ \brainstorm{equation} u(x,0) = I(x), \quad x\in [0,50] \label{wave:pde2:ic:u} \tag{18} \cease{equation} \]
\[ \begin{equation} u_t(10,0) = V(x), \quad ten\in [0,Fifty] \characterization{moving ridge:pde2:ic:ut} \tag{19} \end{equation} \]
\[ \begin{equation} u(0,t) = 0, \quad t>0 \label{moving ridge:pde2:bc:0} \tag{xx} \terminate{equation} \]
\[ \begin{equation} u(Fifty,t) = 0, \quad t>0 \label{wave:pde2:bc:L} \tag{21} \stop{equation} \]
Sampling the PDE at \((x_i,t_n)\) and using the same finite difference approximations every bit above, yields
\[ \begin{equation} [D_tD_t u = c^2 D_xD_x u + f]^{n}_i \label{moving ridge:pde2:fdop} \tag{22} \terminate{equation} \]
Writing this out and solving for the unknown \(u^{northward+1}_i\) results in
\[ \begin{equation} u^{n+1}_i = -u^{due north-ane}_i + 2u^n_i + C^2 (u^{north}_{i+1}-2u^{n}_{i} + u^{due north}_{i-ane}) + \Delta t^2 f^n_i \label{moving ridge:pde2:step3b} \tag{23} \end{equation} \]
The equation for the commencement time step must be rederived. The discretization of the initial condition \(u_t = Five(x)\) at \(t=0\) becomes
\[ [D_{2t}u = V]^0_i\quad\Rightarrow\quad u^{-1}_i = u^{ane}_i - two\Delta t V_i, \]
which, when inserted in (23) for \(due north=0\), gives the special formula
\[ \begin{equation} u^{one}_i = u^0_i + \Delta t V_i + {\frac{i}{2}} C^2 \left(u^{0}_{i+ane}-2u^{0}_{i} + u^{0}_{i-1}\right) + \frac{1}{2}\Delta t^2 f^0_i \characterization{wave:pde2:step3c} \tag{24} \end{equation} \]
Using an belittling solution of concrete significance¶
Many wave problems characteristic sinusoidal oscillations in time and space. For example, the original PDE problem (1)-(5) allows an verbal solution
\[ \brainstorm{equation} u_e(x,t) = A\sin\left(\frac{\pi}{L}x\right) \cos\left(\frac{\pi}{50}ct\right) \label{wave:pde2:test:ue} \tag{25} \terminate{equation} \]
This \(u_e\) fulfills the PDE with \(f=0\), boundary conditions \(u_e(0,t)=u_e(50,t)=0\), every bit well as initial conditions \(I(x)=A\sin\left(\frac{\pi}{L}x\right)\) and \(V=0\).
How to utilise exact solutions for verification.
Information technology is mutual to use such exact solutions of physical interest to verify implementations. Withal, the numerical solution \(u^n_i\) will but be an approximation to \(u_e(x_i,t_n)\). We accept no knowledge of the precise size of the mistake in this approximation, and therefore we can never know if discrepancies between \(u^n_i\) and \(u_e(x_i,t_n)\) are caused by mathematical approximations or programming errors. In item, if plots of the computed solution \(u^n_i\) and the exact one (25) look similar, many are tempted to claim that the implementation works. Nevertheless, even if color plots expect nice and the accuracy is "deemed good", at that place tin can all the same be serious programming errors nowadays!
The only way to apply exact physical solutions like (25) for serious and thorough verification is to run a series of simulations on finer and finer meshes, measure out the integrated mistake in each mesh, and from this data estimate the empirical convergence rate of the method.
An introduction to the calculating of convergence rates is given in Section 3.1.6 in [Langtangen_decay]. At that place is likewise a detailed example on calculating convergence rates in the verification section of the Vibration ODEs affiliate.
In the present problem, ane expects the method to accept a convergence rate of 2 (see the section Analysis of the deviation equations), then if the computed rates are close to two on a sufficiently fine mesh, we take good evidence that the implementation is free of programming mistakes.
Manufactured solution and estimation of convergence rates¶
Specifying the solution and computing corresponding information¶
One problem with the verbal solution (25) is that it requires a simplification (\({V}=0, f=0\)) of the implemented problem (17)-(21). An reward of using a manufactured solution is that we can test all terms in the PDE problem. The idea of this approach is to prepare some chosen solution and fit the source term, boundary conditions, and initial conditions to be uniform with the called solution. Given that our boundary conditions in the implementation are \(u(0,t)=u(L,t)=0\), we must cull a solution that fulfills these conditions. Ane example is
\[ u_e(x,t) = x(50-10)\sin t \]
Inserted in the PDE \(u_{tt}=c^2u_{xx}+f\) we get
\[ -x(L-x)\sin t = -c^2 2\sin t + f\quad\Rightarrow f = (2c^ii - x(L-ten))\sin t \]
The initial weather become
\[\begin{separate} \begin{marshal*} u(ten,0) =& I(x) = 0,\\ u_t(10,0) &= V(x) = ten(L-x) \end{marshal*} \cease{split}\]
Defining a single discretization parameter¶
To verify the code, we compute the convergence rates in a series of simulations, letting each simulation use a finer mesh than the previous one. Such empirical estimation of convergence rates relies on an assumption that some measure out \(E\) of the numerical error is related to the discretization parameters through
\[ E = C_t\Delta t^r + C_x\Delta ten^p, \]
where \(C_t\), \(C_x\), \(r\), and \(p\) are constants. The constants \(r\) and \(p\) are known as the convergence rates in time and space, respectively. From the accuracy in the finite departure approximations, we look \(r=p=two\), since the error terms are of order \(\Delta t^2\) and \(\Delta x^2\). This is confirmed by truncation mistake analysis and other types of analysis.
By using an exact solution of the PDE problem, nosotros will next compute the error measure \(Due east\) on a sequence of refined meshes and meet if the rates \(r=p=two\) are obtained. We volition non exist concerned with estimating the constants \(C_t\) and \(C_x\), merely because we are non interested in their values.
mathcal{I}_t is advantageous to introduce a single discretization parameter \(h=\Delta t=\hat c \Delta 10\) for some constant \(\lid c\). Since \(\Delta t\) and \(\Delta ten\) are related through the Courant number, \(\Delta t = C\Delta ten/c\), nosotros set up \(h=\Delta t\), and and so \(\Delta 10 = hc/C\). Now the expression for the mistake mensurate is greatly simplified:
\[ Due east = C_t\Delta t^r + C_x\Delta ten^r = C_t h^r + C_x\left(\frac{c}{C}\right)^r h^r = Dh^r,\quad D = C_t+C_x\left(\frac{c}{C}\correct)^r \]
Computing errors¶
We choose an initial discretization parameter \(h_0\) and run experiments with decreasing \(h\): \(h_i=2^{-i}h_0\), \(i=i,2,\ldots,m\). Halving \(h\) in each experiment is not necessary, but it is a common choice. For each experiment nosotros must record \(Eastward\) and \(h\). Standard choices of fault measure are the \(\ell^two\) and \(\ell^\infty\) norms of the mistake mesh function \(e^n_i\):
\[ \begin{equation} E = ||e^n_i||_{\ell^2} = \left( \Delta t\Delta x \sum_{n=0}^{N_t}\sum_{i=0}^{N_x} (e^n_i)^2\right)^{\frac{one}{2}},\quad e^n_i = u_e(x_i,t_n)-u^n_i, \label{wave:pde2:fd:MMS:E:l2} \tag{26} \terminate{equation} \]
\[ \begin{equation} E = ||e^n_i||_{\ell^\infty} = \max_{i,n} |east^n_i| \characterization{wave:pde2:fd:MMS:E:linf} \tag{27} \end{equation} \]
In Python, one can compute \(\sum_{i}(e^{n}_i)^two\) at each time step and accumulate the value in some sum variable, say e2_sum . At the final time step one can exercise sqrt(dt*dx*e2_sum) . For the \(\ell^\infty\) norm one must compare the maximum error at a time level ( e.max() ) with the global maximum over the time domain: e_max = max(e_max, due east.max()) .
An alternative mistake measure out is to use a spatial norm at one time step only, east.m., the finish fourth dimension \(T\) (\(northward=N_t\)):
\[ \begin{equation} East = ||due east^n_i||_{\ell^ii} = \left( \Delta x\sum_{i=0}^{N_x} (e^n_i)^2\right)^{\frac{1}{2}},\quad due east^n_i = u_e(x_i,t_n)-u^n_i, \characterization{_auto6} \tag{28} \end{equation} \]
\[ \begin{equation} E = ||e^n_i||_{\ell^\infty} = \max_{0\leq i\leq N_x} |e^{n}_i| \label{_auto7} \tag{29} \terminate{equation} \]
The important point is that the error measure (\(Due east\)) for the simulation is represented by a single number.
Computing rates¶
Allow \(E_i\) exist the error measure in experiment (mesh) number \(i\) (non to be confused with the spatial index \(i\)) and permit \(h_i\) exist the respective discretization parameter (\(h\)). With the error model \(E_i = Dh_i^r\), we can estimate \(r\) by comparison two consecutive experiments:
\[\begin{carve up} \begin{align*} E_{i+1}& =D h_{i+ane}^{r},\\ E_{i}& =D h_{i}^{r} \end{align*} \end{split up}\]
Dividing the 2 equations eliminates the (uninteresting) constant \(D\). Thereafter, solving for \(r\) yields
\[ r = \frac{\ln E_{i+1}/E_{i}}{\ln h_{i+1}/h_{i}} \]
Since \(r\) depends on \(i\), i.e., which simulations we compare, we add an alphabetize to \(r\): \(r_i\), where \(i=0,\ldots,m-two\), if we have \(1000\) experiments: \((h_0,E_0),\ldots,(h_{grand-i}, E_{m-one})\).
In our present discretization of the wave equation nosotros wait \(r=2\), and hence the \(r_i\) values should converge to ii as \(i\) increases.
Constructing an exact solution of the discrete equations¶
With a manufactured or known belittling solution, as outlined higher up, we can approximate convergence rates and encounter if they have the right asymptotic behavior. Feel shows that this is a quite practiced verification technique in that many common bugs will destroy the convergence rates. A significantly better test though, would be to check that the numerical solution is exactly what it should be. This will in general require exact knowledge of the numerical fault, which we do non unremarkably accept (although we in the department Assay of the difference equations establish such noesis in simple cases). Nonetheless, it is possible to look for solutions where we can show that the numerical error vanishes, i.eastward., the solution of the original continuous PDE trouble is as well a solution of the detached equations. This property often arises if the exact solution of the PDE is a lower-guild polynomial. (Truncation fault analysis leads to error measures that involve derivatives of the exact solution. In the present trouble, the truncation error involves 4th-order derivatives of \(u\) in space and time. Choosing \(u\) every bit a polynomial of degree 3 or less will therefore lead to vanishing error.)
We shall now illustrate the construction of an exact solution to both the PDE itself and the discrete equations. Our chosen manufactured solution is quadratic in space and linear in fourth dimension. More specifically, we set
\[ \begin{equation} u_e (x,t) = x(L-x)(1+{\frac{1}{two}}t), \label{moving ridge:pde2:fd:verify:quadratic:uex} \tag{thirty} \terminate{equation} \]
which past insertion in the PDE leads to \(f(x,t)=2(1+t)c^2\). This \(u_e\) fulfills the purlieus weather \(u=0\) and demands \(I(x)=ten(L-x)\) and \(V(x)={\frac{one}{2}}x(50-x)\).
To realize that the called \(u_e\) is also an exact solution of the discrete equations, we first remind ourselves that \(t_n=n\Delta t\) and so that
\[ \begin{equation} \lbrack D_tD_t t^2\rbrack^northward = \frac{t_{n+1}^2 - 2t_n^two + t_{n-1}^2}{\Delta t^two} = (n+1)^2 -2n^ii + (n-1)^2 = two, \label{_auto8} \tag{31} \end{equation} \]
\[ \begin{equation} \lbrack D_tD_t t\rbrack^n = \frac{t_{n+1} - 2t_n + t_{n-ane}}{\Delta t^2} = \frac{((due north+1) -2n + (n-ane))\Delta t}{\Delta t^two} = 0 \characterization{_auto9} \tag{32} \end{equation} \]
Hence,
\[ [D_tD_t u_e]^n_i = x_i(L-x_i)[D_tD_t (1+{\frac{i}{two}}t)]^n = x_i(Fifty-x_i){\frac{ane}{2}}[D_tD_t t]^n = 0 \]
Similarly, we get that
\[\begin{split} \begin{align*} \lbrack D_xD_x u_e\rbrack^n_i &= (1+{\frac{1}{ii}}t_n)\lbrack D_xD_x (xL-x^2)\rbrack_i\\ & = (1+{\frac{one}{2}}t_n)\lbrack LD_xD_x x - D_xD_x x^2\rbrack_i \\ &= -ii(1+{\frac{i}{2}}t_n) \end{align*} \stop{split}\]
At present, \(f^n_i = two(1+{\frac{i}{2}}t_n)c^two\), which results in
\[ [D_tD_t u_e - c^2D_xD_xu_e - f]^n_i = 0 + c^2 2(ane + {\frac{1}{two}}t_{northward}) + 2(1+{\frac{1}{2}}t_n)c^2 = 0 \]
Moreover, \(u_e(x_i,0)=I(x_i)\), \(\fractional u_e/\partial t = V(x_i)\) at \(t=0\), and \(u_e(x_0,t)=u_e(x_{N_x},0)=0\). Also the modified scheme for the first fourth dimension footstep is fulfilled by \(u_e(x_i,t_n)\).
Therefore, the exact solution \(u_e(ten,t)=x(L-x)(i+t/2)\) of the PDE problem is likewise an exact solution of the detached problem. This means that we know beforehand what numbers the numerical algorithm should produce. Nosotros tin can utilize this fact to check that the computed \(u^n_i\) values from an implementation equals \(u_e(x_i,t_n)\), within machine precision. This result is valid regardless of the mesh spacings \(\Delta x\) and \(\Delta t\)! Nevertheless, there might exist stability restrictions on \(\Delta ten\) and \(\Delta t\), so the test can merely be run for a mesh that is compatible with the stability benchmark (which in the present example is \(C\leq 1\), to exist derived subsequently).
Notice.
A product of quadratic or linear expressions in the various independent variables, as shown above, will often fulfill both the PDE problem and the discrete equations, and can therefore be very useful solutions for verifying implementations.
Withal, for 1D wave equations of the blazon \(u_{tt}=c^2u_{xx}\) we shall see that there is always some other much more powerful way of generating verbal solutions (which consists in just setting \(C=one\) (!), as shown in the section Analysis of the divergence equations).
Waves On A String Simulation,
Source: https://www.devitoproject.org/devito_book/notebooks/02_wave/wave1D_fd1.html
Posted by: garnerlifivend1962.blogspot.com
0 Response to "Waves On A String Simulation"
Post a Comment