Using finite difference to solve DE

1. Jul 29, 2016

jdawg

1. The problem statement, all variables and given/known data
d2T/dx2 = 5*(dT/dx) - 0.1*x = 0
T(0) = 50
T(10) = 400
(Δx) = 2

I've figured out how to do these problems when Δx = 1, but when it equals any other number it goes wrong.
I know you start by plugging in the algebraic approximations for the differential elements, I think maybe my problem is the nodes?

2. Relevant equations
d2T/dx2 = (Ti+1 * Ti + Ti-1)/((Δx))
dT/dx = (Ti+1 - Ti-1)/(2*Δx)

3. The attempt at a solution
node1 = 2; node2 = 4; node3 = 6; node4 = 8;

For node1:
(T2 - 2*T1 + T0)/(22) + 5*((T2 - T0)/(2*2)) - 0.1*(2) = 0

End up with ======> -0.5*T1 + 1.5*T2 = 50.2

For node2:
(T3 - 2*T2 + T1)/(22) + 5*((T3 - T1)/(2*2)) - 0.1*(4) = 0

End up with =======> -T1 - 0.25*T2 + 1.5*T3 = 0.4

For node3:
(T4 - 2*T3 + T2)/(22) + 5*((T4 - T2)/(2*2)) - 0.1*(6) = 0

End up with ========> -T2 - 0.5*T3 + 1.5*T4 = 0.6

For node4:
(T5 - 2*T4 + T3)/(22) + 5*((T5 - T3)/(2*2)) - 0.1*(8) = 0

End up with ========> -T3 - 0.5*T4 = -724.2

I put all the coefficients into a matrix and its tridiagonal, which is good I think. But when I try to plot it in matlab it gives me a crazy looking zigzag, which I'm pretty sure isn't correct.

If someone could point out what I'm doing wrong I would really appreciate it! :)

Last edited by a moderator: Jul 30, 2016
2. Jul 30, 2016

Simon Bridge

For arbitrary $\Delta x$ you have: $x_n = x_0 + n\Delta x$ ... this is the lynchpin.
Tri-diagonal matrix is correct: that's what you get for a second derivative.

If you only had the second derivative, then the matrix would have (1,-2,1) on the diagonal, multiplied by a constant.
You should be able to write out the matrices for each term and just add them up.

3. Jul 30, 2016

jdawg

I still don't really understand.

4. Jul 31, 2016

Simon Bridge

It is best to understand the techniques you are applying instead of just going through the motions by rote.

Note: post #1, under "relevant equations", check that first expression.

If $D^2_iy$ indicates the second derivative of $y(x)$ evaluated at point $x=x_i : i\in\mathbb{Z}$ then:
$$D^2_iy \approx \frac{y_{i+1}-2y_{i}+y_{i-1}}{h^2}: y_i=y(x_i)$$
... here I used $h=\Delta x$ to avoid typing lots of deltas.

However, this expression may be inconsistent with the second expression in the same section.
If $D_iy\approx\frac{1}{2h}\big(y_{i+1}-y_{i-1}\big)$ then: $D_i^2y \approx \frac{1}{4h}\big(D_{i+1}y-D_{i-1}y\big)$
... so it kinda helps to know where you are coming from.

This would make the matrix for a second order DE potentially quin-diagonal. But I don't think it's a fatal problem since these are only approximations anyway.

Now, using this notation...
I can define a vector $\vec x = \big(x_0,x_1,x_2,\cdots x_n,\cdots x_{N-2},x_{N-1}\big)^t : x_n=x_0+nh$ ... for the N discrete points you want to do your evaluation over. (The "$^t$" indicates a transpose so this is defined to be a column vector.)

I can derive, from that, other vectors like: $\vec y = \big(y_0,y_1, \cdots y_n, cdots y_{N-1}\big)^t:y_n=y(x_n)$ and, similarly $\vec y' = \big(y'(x_0)\cdots\big)^t$ etc.

When I look at the relationships, I can see that $\vec y' = D\vec y$ where $D$ is a matrix that looks like this:
$$D=\frac{1}{2h}\begin{bmatrix} 0 & 1 & 0 &\cdots\\ -1 & 0 & 1 & 0 & \cdots\\ 0 & -1 & 0 & 1 & 0 & \cdots\\ \cdots & 0 & -1 & 0 & 1 & 0 & \cdots\\ & & \ddots & \ddots & \ddots & \ddots & \ddots & &\\ \end{bmatrix}$$ ... ie, this is a tri-diagonal matrix where the diagonal elements are (-1,0,1) times a constant (1/2h).

This follows from your definition in post #1 (section 2, 2nd equation). You should check, if you have not already.
OTOH: If we follow the definition of the derivative given in math text books you get a tri-diagonal matrix with elements (0,-1,1)/h.

... however you do it, it also follows that $\vec y'' = D\vec y' = D(D\vec y) = D^2\vec y$

Thus any DE can, with a bit of care, be rewritten as a matrix equation.

Yours is (from post #1):
Which is $D^2\vec T = 5D\vec T - 0.1\vec x$
Where $h=2, T(0)=T(x_0)=T_0=50, T(10)=T(x_{175})=T_{175} = 400$

Rearrange, you get $A\vec T = \vec x$ where $A=50D-10D^2$.
The solution, therefore, amounts to finding the inverse of the matrix A.

Last edited: Jul 31, 2016