CODEWITHSUNDEEP
algorithmnumerical-methodsodescientific-computing
Insight
Insight

Reputation: 31

Finite Difference Method for Solving ODEs Algorithm

I'm trying to devise an algorithm for the finite difference method, but I'm a bit confused. The ODE in question is y''-5y'+10y = 10x, with y(0)=0 and y(1)=100. So I need a way to somehow obtain the coefficients that will multiply "y_i" from the relation:

enter image description here

And then store the resultant coefficients into a matrix, which will be the matrix of the system I'll solve trough Gauss-Jordan. The question boils down to how to obtain these coefficients and move them to the matrix. I thought about working out the coefficients by hand and then just inputing the matrix, but I need to do this for steps of size 0.1, 0.001 and 0.001, so that's really not a viable option here.

Upvotes: 0

Views: 190

Answers (1)

Nico Schertler
Nico Schertler

Reputation: 32587

Let us assume the more general case of the ODE

c1 * y''(x) + c2 * y'(x) + c3 * y(x) + c4 * x = 0

with the boundary conditions

y(0) = lo
y(1) = hi

And you want to solve this for x ∈ [0, 1] with a step size of h = 1 / n (where n + 1 is the number of samples). We want to solve for the yi = y(h * i). The yi range from i ∈ [0, n]. To do this, we want to solve a linear system

A y = b

Every interior yi will impose a single linear constraint. Hence, we have n - 1 rows in A and n - 1 columns corresponding to the unknown yi.

To set up A and b, we can simply slide a window over our unknown yi (I assume zero-based indexing).

A = 0  //the zero matrix
b = 0  //the zero vector
for i from 1 to n - 1
    //we are going to create the constraint for yi and store it in row i-1

    //coefficient for yi+1
    coeff = c1 / h^2 + c2 / h
    if i + 1 < n
        A(i - 1, i) = coeff
    else
        b(i - 1) -= coeff * hi  //we already know yi+1

    //coefficient for yi
    coeff = -2 * c1 / h^2 - c2 / h + c3
    A(i - 1, i - 1) = coeff

    //coefficient for yi-1
    coeff = c1 / h^2
    if i - 1 > 0
        A(i - 1, i - 2) = coeff
    else
        b(i - 1) -= coeff * lo  //we already know yi-1

    //coefficient for x
    b(i - 1) -= c4 * i * h
 next

Upvotes: 2

Related Questions