Using finite difference to solve DE

[jdawg]

Homework Statement

d²T/dx² = 5*(dT/dx) - 0.1*x = 0
T(0) = 50
T(10) = 400
(Δx) = 2I'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?

Homework Equations

d²T/dx² = (T_i+1 * T_i + T_i-1)/(^(Δx))
dT/dx = (T_i+1 - T_i-1)/(2*Δx)

The Attempt at a Solution

node1 = 2; node2 = 4; node3 = 6; node4 = 8;

For node1:
(T₂ - 2*T₁ + T₀)/(2²) + 5*((T₂ - T₀)/(2*2)) - 0.1*(2) = 0

End up with ======> -0.5*T₁ + 1.5*T₂ = 50.2

For node2:
(T₃ - 2*T₂ + T₁)/(2²) + 5*((T₃ - T₁)/(2*2)) - 0.1*(4) = 0

End up with =======> -T₁ - 0.25*T₂ + 1.5*T₃ = 0.4

For node3:
(T₄ - 2*T₃ + T₂)/(2²) + 5*((T₄ - T₂)/(2*2)) - 0.1*(6) = 0

End up with ========> -T₂ - 0.5*T₃ + 1.5*T₄ = 0.6

For node4:
(T₅ - 2*T₄ + T₃)/(2²) + 5*((T₅ - T₃)/(2*2)) - 0.1*(8) = 0

End up with ========> -T₃ - 0.5*T₄ = -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! :)

 

[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.

 

[jdawg]

I still don't really understand.

 

[Simon Bridge]

I still don't really understand.

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.

You follow me so far?

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 textbooks 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):

d2T/dx2 = 5*(dT/dx) - 0.1*x = 0

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.