Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Second order system of odes with variable coefficients

  1. Apr 28, 2014 #1
    Hi,

    I have looked everywhere. Can someone please point me in the right direction for solving a system of ODEs with variable coefficients? I managed to solve such system with constant coefficients.
     
  2. jcsd
  3. Apr 28, 2014 #2

    lurflurf

    User Avatar
    Homework Helper

    Solving symbolically or numerically?
    Symbolically there is the Peano-Baker Method and use of integrating factors.
    it helps if the functions are simple
    Have functions in mind?
    numerically variable coefficients are no problem if they are well behaved

    Have you worked out simple examples like
    $$\begin{cases} \mathrm{x}^{\prime}(t)=-t\,\mathrm{y}(t) \\ \mathrm{y}^{\prime}(t)=\phantom{-}t\,\mathrm{x}(t) \end{cases}$$
    or
    $$\begin{cases} \mathrm{x}^{\prime}(t)=t\,\mathrm{y}(t) \\ \mathrm{y}^{\prime}(t)=\phantom{t}\,\mathrm{x}(t) \end{cases}$$

    you can see how quickly things get messy
    http://journals.cambridge.org/downl...52a.pdf&code=5491acce57881b30db323feb4fadd1e3
     
  4. May 1, 2014 #3
    Thank you. No I haven't worked out those examples. Do you know of any good examples of ODEs that are solved using the Peano Baker method?
     
  5. May 1, 2014 #4
    Thanks, my equations in the A matrix are:

    [ alpha / (b + 2*theta*x) -alpha / (b + 2*theta*x) ]
    [ -alpha / (b + 2*theta*(L-x)) alpha / (b + 2*theta*(L-x)) ]

    where x is the only variable and this gets complicated enough.
     
  6. May 3, 2014 #5
    After several attempts, I am stuck on the second iteration. I tried integrating in Matlab and an integral does not exist, so I am giving up on this one. Do you think power series would work here e.g. Frobenius method?
     
  7. May 4, 2014 #6

    lurflurf

    User Avatar
    Homework Helper

    Do you mean
    $$
    \begin{bmatrix} \\ \mathrm{u}^{\prime}(x) \\ \\ \mathrm{v}^{\prime}(x) \\ \\ \end{bmatrix}
    = \alpha\begin{bmatrix}
    \phantom{-} \dfrac{1}{b+2\theta x} & -\dfrac{1}{b+2\theta x}
    \\
    -\dfrac{1}{b+2\theta (L-x)} & \phantom{-}\dfrac{1}{b+2\theta (L-x)}
    \end{bmatrix}\begin{bmatrix} \\ \mathrm{u}(x) \\ \\ \mathrm{v}(x) \\ \\ \end{bmatrix}$$

    That is an interesting equation
    if it is the one you mean it is easily solved
    it might help to rename some parameters to make the algebra less messy
    write a differential equation for u-v (no u+v) and solve it
    write a differential equation for u+v (u-v is now known) and solve it
    the final answer has an integral of a radical function that can be left unevaluated or evaluated with hypergeometric function

    http://www.wolframalpha.com/input/?.../(b+2theta*(l-t))+alpha*z(t)/(b+2theta*(l-t))
     
    Last edited: May 4, 2014
  8. May 4, 2014 #7

    lurflurf

    User Avatar
    Homework Helper

    $$\begin{cases}\frac{\big(\mathrm{u}(x)-\mathrm{v}(x)\big)^\prime}{\mathrm{u}(x)-\mathrm{v}(x)}=\alpha\left(\frac{1}{b+2\, \theta \, x}+\frac{1}{b+2\, \theta \, (L-x)}\right) \\ \frac{\big(\mathrm{u}(x)+\mathrm{v}(x)\big)^\prime}{\mathrm{u}(x)-\mathrm{v}(x)}=\alpha\left(\frac{1}{b+2\, \theta \, x}-\frac{1}{b+2\, \theta \, (L-x)}\right)\end{cases}$$

    Hopefully no typos
     
  9. May 4, 2014 #8

    lurflurf

    User Avatar
    Homework Helper

    I do not have a textbook handy.
    Besides the above ones some simple examples are
    $$\begin{bmatrix} \mathrm{u}^{\prime}(t) \\ \mathrm{v}^{\prime}(t)\end{bmatrix}=
    \begin{bmatrix} 0 & 1 \\ a & -a \, t\end{bmatrix}\begin{bmatrix} \mathrm{u}(t) \\ \mathrm{v}(t)\end{bmatrix}$$

    $$\begin{bmatrix} \mathrm{u}^{\prime}(t) \\ \mathrm{v}^{\prime}(t)\end{bmatrix}=
    \begin{bmatrix} a & e^{-t} \\ 0 & a \end{bmatrix}\begin{bmatrix} \mathrm{u}(t) \\ \mathrm{v}(t)\end{bmatrix}$$

    http://arxiv.org/pdf/1011.1775v2.pdf
    has two examples

    $$\begin{bmatrix} \mathrm{u}^{\prime}(t) \\ \mathrm{v}^{\prime}(t)\end{bmatrix}=
    \begin{bmatrix} 0 & t \\ a & 0\end{bmatrix}\begin{bmatrix} \mathrm{u}(t) \\ \mathrm{v}(t)\end{bmatrix}$$

    $$\begin{bmatrix} \mathrm{u}^{\prime}(t) \\ \mathrm{v}^{\prime}(t)\end{bmatrix}=
    \begin{bmatrix} 1 & t \\ 0 & a \, t\end{bmatrix}\begin{bmatrix} \mathrm{u}(t) \\ \mathrm{v}(t)\end{bmatrix}$$
     
  10. May 4, 2014 #9
    Thanks.

    Here is my issue, I even tried solving the equation numerically and the solution blows up incrementally using central differencing.

    Here is a simple example system:

    d^2y1/dx^2 = alpha * ( y1 - y2 )
    d^2y2/dx^2 = alpha * ( y2 - y1 )

    As I don't have any initial condition for y(x-h) i.e. y-1 where h is the timestep size, I am trying to employ the second order forward difference method, but this requires two unknowns y(x+2h) and y(x+h). Any idea what to do?
     
  11. May 4, 2014 #10

    lurflurf

    User Avatar
    Homework Helper

    What are the initial conditions?
    Normally one estimates
    y(x+h) and y(x+2h) from the derivative

    This one is another where sum difference is helpful

    (y1-y2)''=2*alpha * ( y1 - y2 )
    (y1+y2)''=0

    The solution has e^x so if x is large unless the coefficient is small or zero the solution tends to rapid growth.
     
  12. May 4, 2014 #11
    ICs:

    y1(0) = 0
    y2(0) = 1000
    I tried alpha small 0.0005 and it still grows rapidly, except that this is not the case with the analytical solution
     
  13. May 4, 2014 #12

    lurflurf

    User Avatar
    Homework Helper

    you need four conditions such as
    y1(0), y2(0), y1'(0),and y2'(0)

    Probably your numerical method is unstable
     
  14. May 5, 2014 #13
    Of course, I should have included the other 2 for you. y1(0)=0, y2(0)=1000, y1(n)=1000, y2(n)=0. However, this does not seem to work with central differencing no matter what I do as a first derivative BC is required to solve for the first y at -1
     
  15. May 5, 2014 #14

    lurflurf

    User Avatar
    Homework Helper

    What is n? So this is a boundary value problem so it is a bit different. Some of the simple options include
    1)Shooting method
    guess y1(1) and y2(1) then you have an initial value problem
    compute y1(n) and y2(n) and compare to the given values
    repeat until y1(n) and y2(n) match the given values to within tolerance
    2)Direct sparse
    write 2n-2 equations in 2n-2 variables then direct solve them
    2)Iterative sparse
    write 2n-2 equations in 2n-2 variables then iterate them
     
  16. May 5, 2014 #15
    To restate the system:

    d^2y1/dx^2 = alpha * ( y1 - y2 )
    d^2y2/dx^2 = alpha * ( y2 - y1 )

    Where alpha is a constant

    In numerical form:

    [ y1(I+1) - 2y1(I) + y1(I-1) ] / h^2 = alpha * [ (y1(I) - y2(I) ]
    [ y2(I+1) - 2y2(I) + y2(I-1) ] / h^2 = alpha * [ (y2(I) - y1(I) ]

    Where h is the timestep size

    Hence,

    y1(I+1) = (2 + alpha * h^2) * y1(I) - y1(I-1) - alpha * h^2 * y2(I)
    y2(I+1) = (2 + alpha * h^2) * y2(I) - y2(I-1) - alpha * h^2 * y1(I)

    Say I want to use the shooting method, with my BCs y1(0)=0 and y2(0)=sigma0 where sigma0 is any constant. So I need an initial y1'(0) and y2'(0). Given the central difference method [ y1(I+1) - y1(I-1) ] / (2h) where h is the timestep. Thus initially, I have:

    a) [ y1(1) - y1(-1) ] / 2h = a1 for my initial guess for the derivative, hence y1(-1) = y1(1) - 2*h*a1
    b) Similarly, y2(-1) = y2(1) - 2*h*a2
    c) When I=0, substituting from (a) and (b) into the above numerical expressions:

    y1(1) = (2 + alpha * h^2) * y1(0) - y1(-1) - alpha * h^2 * y2(0)
    y1(1) + y1(1) - 2*h*a1 = (2 + alpha * h^2) * y1(0) - alpha * h^2 * y2(0)
    y1(1) = [ (2 + alpha * h^2) * y1(0) - alpha * h^2 * y2(0) - 2 * h * a1 ] / 2

    y2(1) = (2 + alpha * h^2) * y2(0) - y2(-1) - alpha * h^2 * y1(0)
    y2(1) + y2(1) - 2*h*a2 = (2 + alpha * h^2) * y2(0) - alpha * h^2 * y1(0)
    y2(1) = [ (2 + alpha * h^2) * y2(0) - alpha * h^2 * y1(0) - 2 * h * a2 ] / 2

    For each timestep at I=1 to N

    % Iterate here until difference between IVP and BVP is small

    y1(I+1) = (2 + alpha * h^2) * y1(I) - y1(I-1) - alpha * h^2 * y2(I)
    y2(I+1) = (2 + alpha * h^2) * y2(I) - y2(I-1) - alpha * h^2 * y1(I)

    End for
     
  17. May 6, 2014 #16
    One of the terms in the analytic solution is exponentially growing. This term could be giving you problems. In order to determine if this is the problem we need to know both the length of the system and your mesh spacing. You want the mesh spacing to be smaller than the characteristic distance of the exponential term. You also don't want the domain size to be much larger than this characteristic distance.

    Also what values are you using for as your initial guesses for [itex]y_1'[/itex] and [itex]y_2'[/itex]?

    It may be worth wild to further simplify the system.
    What happens if you set [itex]\alpha=0[/itex]. Is you system stable? Do reproduce the analytic solution for this case?

    Second, what happens if you try solving the equation
    [itex]y_1''=\alpha y_1[/itex]?
     
    Last edited: May 6, 2014
  18. May 8, 2014 #17
    Thanks for the reply. Does reduction of order work on systems of ODEs with variable coefficients?
     
  19. May 10, 2014 #18

    lurflurf

    User Avatar
    Homework Helper

    ^Yes reduction of order works on systems of ODEs with variable coefficients.
     
  20. May 12, 2014 #19
    I thought so. Here is what I did and the solution didn't pan out.


    d^2(sigma1)/dx^2 = -(2*theta) / (b + 2*theta*x) * d(sigma1)/dx + alpha / (b + 2*theta*x) * (sigma1 - sigma2)
    d^2(sigma2)/dx^2 = (2*theta) / [ b + 2*theta*(L-x) ] * d(sigma2)/dx + alpha / [ b + 2*theta*(L-x) ] * (sigma2 - sigma1)


    1) Reducing these to a first order system:



    Let y1 = sigma1, y2 = d(sigma1)/dx

    d(y1)/dx = d(digma1)/dx = y2 and;

    d(y2)/dx = d^2(sigma1)/dx^2 = -(2*theta) / (b + 2*theta*x) * d(sigma1)/dx + alpha / (b + 2*theta*x) * (sigma1 - sigma2)



    Let y3 = sigma2, y4 = d(sigma2)/dx

    d(y3)/dx = d(sigma2)/dx = y4 and;

    d(y4)/dx = d^2(sigma2)/dx^2 = (2*theta) / [ b + 2*theta*(L-x) ] * d(sigma2)/dx + alpha / [ b + 2*theta*(L-x) ] * (sigma2 - sigma1)



    Then the first order system is:



    d(y1)/dx = y2

    d(y2)/dx = -(2*theta) / (b + 2*theta*x) * y2 + alpha / (b + 2*theta*x) * (y1 - y3)

    d(y3)/dx = y4

    d(y4)/dx = (2*theta) / [ b + 2*theta*(L-x) ] * y4 + alpha / [ b + 2*theta*(L-x) ] * (y3 - y1)



    I then take the Laplace transforms, applying the BCs (y1(0) = sigma1(0) = 0, y3(0) = sigma2(0) = sigma0) to get:



    s*Y1(s) - Y2(s) = 0

    -alpha*Y1(s) + b*s*Y2(s) + alpha*Y3(s) = b*y2(0)

    s*Y3(s) - Y4(s) = sigma0

    alpha*Y1(s) - alpha*Y3(s) + (bs + 2*theta*L*s)*Y4(s) = (b + 2*theta*L)*y4(0)



    Then I solve the system using Gauss elimination and back substitute to make sure I got the correct answer, which I did.



    Finally, I take the inverse Laplace transform in Matlab and got a solution for y1 and y3 which represent sigma1 and sigma2 respectively. Differentiating them appropriately and plugging them back into the original ODEs does not get me the correct answer.



    Even though my boundary conditions for the first derivative are unknown at x=0, I would still expect that these steps yield a correct solution because the unknowns at x=0 are still constants y2(0) and y4(0).



    I would appreciate any suggestions.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Second order system of odes with variable coefficients
Loading...