Scribd
Upload
14 views7 pages

Finite Difference Method for ODEs

finite difference method- matrix method solution for Machine Learning

Uploaded by

Disha Sen
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
14 views7 pages

Finite Difference Method for ODEs

finite difference method- matrix method solution for Machine Learning

Uploaded by

Disha Sen
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

Finite Difference Method:

Another way to solve the ODE boundary value problems is the finite difference method, where
we can use finite difference formulas at evenly spaced grid points to approximate the differential
equations. This way, we can transform a differential equation into a system of algebraic
equations to solve.
In the finite difference method, the derivatives in the differential equation are approximated
using the finite difference formulas. We can divide the interval of [a,b] into n equal subintervals
of length h as shown in the following figure.

Commonly, we usually use the central difference formulas in the finite difference methods due to
the fact that they yield better accuracy. The differential equation is enforced only at the grid
points, and the first and second derivatives are:

𝑑𝑦 𝑦𝑖+1 − 𝑦𝑖−1
=
𝑑𝑥 2ℎ
𝑑2 𝑦 𝑦𝑖−1 − 2𝑦𝑖 + 𝑦𝑖+1
=
𝑑𝑥 2 ℎ2
These finite difference expressions are used to replace the derivatives of y in the differential
equation which leads to a system of n+1 linear algebraic equations if the differential equation is
linear. If the differential equation is nonlinear, the algebraic equations will also be nonlinear.

Example: Consider the following second order equation of rocket motion:


𝑑2𝑦
= −𝑔
𝑑𝑥 2
Boundary conditions are: y(0)=0, y(10)=50. Let’s take n=50. Since the time interval is [0,
10] and we have n=50, therefore, h=0.2, using the finite difference approximated derivatives, we
have:
𝑦0 = 0
𝑦𝑖−1 − 2𝑦𝑖 + 𝑦𝑖+1 = −𝑔ℎ2 , 𝑖 = 1, 2, 3, … … … , 𝑛 − 1
𝑦10 = 50
In matrix form: AY=B, the components of A are basically the co-efficient of yi for each equations.
For first row put i=0, then coefficients of y0 is 1 and other terms of this row (here, rest of the 49
terms) are all 0. Similarly it will happen for other equations also.

Therefore, we have n+1=51 equations in the system.


Python Script:
Eigen Value Problem:
Example: Consider the 1D Schrödinger equation for a particle in a box (potential energy,
V=0)
𝑑2𝜓 2𝑚
= − 𝐸𝜓
𝑑𝑥 2 ℏ2
We discretize the domain (width of the potential) into equal step size: x=h such that we can
write:
𝜓𝑖−1 −2𝜓𝑖 +𝜓𝑖+1 2𝑚
=− 𝐸𝜓𝑖
ℎ2 ℏ2

h is just step size, don’t mix with Planck’s constant.


𝜓0 −2𝜓1 +𝜓2 2𝑚
Put i=1, =− 𝐸𝜓1
ℎ2 ℏ2
𝜓1 −2𝜓2 +𝜓3 2𝑚
Put i=2, =− 𝐸𝜓1
ℎ2 ℏ2

……………………………
𝜓𝑛−1 −2𝜓𝑛 +𝜓𝑛+1 2𝑚
Put i=n, =− 𝐸𝜓1
ℎ2 ℏ2

In matrix form:
Let’s define the grid points xi such that:
𝑥1 = 0, 𝑥2 = ℎ, ……, 𝑥𝑛 = (𝑛 − 1)ℎ with boundary conditions 𝜓(0)= 𝜓(𝐿)=0
We will construct a matrix A such that:
2𝑚𝐸
𝐴𝜓 = 𝜓
ℏ2
This matrix A represents the finite difference approximation of the second derivative. This will
be a tridiagonal matrix, that means except the diagonal, top and bottom diagonal elements are
zero. The diagonal elements represent −2⁄ℎ2 , and the off-diagonal elements represent 1⁄ℎ2 .
The matrix looks like this for N grid points:
Therefore,

= E

Let’s understand through python script step by step:

Here we are basically discretizing the x using linspace. If we put just ‘N’ in linspace then the
division is little bit problem as shown in above. So it’s better to take ‘N+1’ in linspace as
following:

Now let’s form the matrix A:


But, here my diagonal element for matrix A is -2. Then we have to just multiply -2 with [Link].

Next, to make the elements for upper diagonal ‘[Link]([Link](N-2),1)’ and lower diagonal
‘[Link]([Link](N-2),-1)’
To check the dimension of the matrix:

ℏ2 𝐸
Now, we will multiply the constant with matrix A and form a new matrix H:
2𝑚ℎ2

So, we have now 𝐻𝜓 = 𝐸𝜓


We want to solve the eigenvalue problem. The inbuilt function of numpy will give us directly the
eigenvalues i.e energies, here and the eigenvectors i.e. 𝜓s.
Thus, we also can plot probability density by plotting 𝜓 ∗ 𝜓.

To solve the s-wave Schrödinger equation for the ground and first excited states of the
hydrogen atom using the matrix method, we first need to discretize the problem and then solve it
numerically using matrix diagonalization.
Schrödinger Equation for Hydrogen Atom (s-wave):
The radial Schrödinger equation for the hydrogen atom (with l=0l = 0l=0, i.e., s-wave) is given
by:
ℏ2 𝑑 2 𝑒 2
[− − ] 𝑢(𝑟) = 𝐸𝑢(𝑟)2
2𝑚 𝑑𝑟 2 𝑟

Where, 𝑢(𝑟) is the radial wavefunction


𝑢(𝑟)

Common questions

Powered by AI

The matrix A is constructed to solve the eigenvalue problem in the Schrödinger equation by using a finite difference approximation of the second derivative. The diagonal elements of matrix A are set to −2/ℎ2, while the off-diagonal elements are set to 1/ℎ2. For N grid points, this forms a tridiagonal matrix where only the diagonal, top, and bottom-diagonal elements are non-zero. The matrix A is then multiplied by the factor ℏ2𝐸/2𝑚ℎ2 to form the matrix H such that Hψ = Eψ, where ψ represents the eigenvector solutions .

The tridiagonal matrix in finite difference approximations is constructed with specific elements: diagonal elements are set to −2/ℎ2, and off-diagonal elements to 1/ℎ2, forming a band structure with zero components elsewhere. This configuration efficiently represents second derivative discretizations in numerical methods (like for Schrödinger or other differential equations), significantly reducing computational effort while maintaining accuracy. The tridiagonal form ensures that the resulting system of algebraic equations is sparse and can be solved efficiently using tailored algorithms .

In the rocket motion example, the following initial setup is given: the second-order differential equation 𝑑2𝑦/𝑑𝑥2 = −𝑔, with boundary conditions y(0)=0 and y(10)=50. The interval is divided into 50 equal subintervals (n=50), resulting in a step size of h=0.2. Using the finite difference method, the differential equation is transformed into a system of algebraic equations with n+1=51 equations to solve .

For the second-order differential equation of rocket motion, 𝑑2𝑦/𝑑𝑥2 = −𝑔, the boundary conditions are y(0)=0 and y(10)=50. The interval [0, 10] is divided into 50 equal subintervals (n=50) with a step size h=0.2. Applying the finite difference approximation to 𝑑2𝑦/𝑑𝑥2 results in 𝑦𝑖−1 −2𝑦𝑖 + 𝑦𝑖+1 = −𝑔ℎ2 for i=1, 2, ..., n−1. The first and last equations represent the boundary conditions: y₀ = 0 and y₁₀ = 50. This forms a system of 51 linear algebraic equations that can be represented in matrix form as AY = B, where A includes the coefficients of yi for the equations .

Central difference formulas are often favored in finite difference methods due to their higher accuracy compared to forward or backward difference schemes. They approximate the derivatives more accurately by considering values at both sides of the grid point where the derivative is being evaluated, minimizing error and providing better representation of the differential equation at those points .

The finite difference method approximates differential equations by using finite difference formulas to replace the derivatives at evenly spaced grid points. This transforms the differential equation into a system of algebraic equations. For instance, using central difference formulas provides approximations for the first and second derivatives: 𝑑𝑦/𝑑𝑥 ≈ (𝑦𝑖+1 − 𝑦𝑖−1)/(2ℎ) and 𝑑2𝑦/𝑑𝑥2 ≈ (𝑦𝑖−1 − 2𝑦𝑖 + 𝑦𝑖+1)/ℎ2. These expressions replace the derivatives in the differential equation, leading to a system of linear algebraic equations for linear differential equations, or nonlinear algebraic equations for nonlinear differential equations .

Matrix diagonalization plays a crucial role in solving the Schrödinger equation for the hydrogen atom by transforming the problem into an eigenvalue problem. Once the matrix H, representing the finite difference approximation of the equation, is set up, diagonalization enables the determination of eigenvalues which correspond to energy levels (ground and first excited states) and eigenvectors which provide the corresponding wavefunctions. This numerical approach simplifies the computation and directly provides the required quantum mechanical solutions .

In the finite difference method, when dealing with nonlinear differential equations, the resulting algebraic equations will also be nonlinear as opposed to linear equations derived from linear differential equations. This is because the finite difference approximations used to replace the derivatives inherently preserve the nature of the differential equations they approximate, requiring nonlinear solution strategies or iterative methods to solve the system of nonlinear equations .

In the matrix form solution for the rocket motion differential equation example, boundary conditions are integrated directly into the matrix formulation. For instance, for y(0) = 0, the first row of matrix A in the system of equations will contain 1 for the coefficient of y₀ and 0 for all other coefficients, effectively setting y₀ equal to 0. Similarly, for y(10) = 50, the last equation directly reflects this condition. These constraints ensure that the solution respects the specified boundary values while solving the matrix system AY=B .

Using 'N+1' in linspace is recommended when discretizing the domain to ensure that the division of the interval includes both endpoints, making it more accurate and evenly distributed. This approach captures the boundary conditions effectively, which is essential for setting up the matrix form solution correctly in the Schrödinger equation .

You might also like