1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Using finite difference to solve DE

  1. Jul 29, 2016 #1
    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. jcsd
  3. Jul 30, 2016 #2

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    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.
     
  4. Jul 30, 2016 #3
    I still don't really understand.
     
  5. Jul 31, 2016 #4

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    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 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
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted