I'm solving a system of stiff ODEs, at first I wanted to implement BDF, but it seem to be a quite complicated method, so I decided to start with Backward Euler method. Basically it says that you can solve an ODE: y ′ = f ( t , y ), with initial guess y ( t 0 ) = y 0 Using the following approximation: y k + 1 = y k + h f ( t k + 1 , y k + 1 ), where h is a step size on parameter tWikipedia article says that you can solve this equation using Newton-Raphson method, which is basically a following iteration: x n + 1 = x n − g ( x n ) g ′ ( x n ) So, the question is how to correctly mix them together? What initial guess x 0 and function g should be?Also f is quite complex in my case and I'm not sure if it possible to find another derivative of it analytically. I want to write an implementation of it by myself, so pblackefined Mathematica functions wouldn't work.

Question

I&#039;m solving a system of stiff ODEs, at first I wanted to implement BDF, but it seem to be a quite complicated method, so I decided to start with Backward Euler method. Basically it says that you can solve an ODE:      y    ′    =  f  (  t  ,  y  ), with initial guess   y  (      t    0    )  =      y    0  Using the following approximation:      y          k      +      1        =      y    k    +  h  f  (      t          k      +      1        ,      y          k      +      1        ), where   h is a step size on parameter   tWikipedia article says that you can solve this equation using Newton-Raphson method, which is basically a following iteration:      x          n      +      1        =      x    n    −            g      (              x        n            )                      g        ′            (              x        n            )      So, the question is how to correctly mix them together? What initial guess       x    0   and function   g should be?Also   f is quite complex in my case and I&#039;m not sure if it possible to find another derivative of it analytically. I want to write an implementation of it by myself, so pblackefined Mathematica functions wouldn&#039;t work.

Nancy Phillips · Accepted Answer

Your aim is to solve the following equation for       y          k      +      1      :  g  (      y          k      +      1        )  :=      y          k      +      1        −      y    k    −  h  f  (      t          k      +      1        ,      y          k      +      1        )  =  0  ,where   f is a known function and       y    k    ,      t          k      +      1       and   h are known values. This gives you the   g for a Newton-Raphson method. As an initial guess, I&#039;d suggest that you use a 1-step forward Euler method to explicitly calculate                    y        ^                    k      +      1        :=      y    k    +  h  f  (      t          k      +      1        ,      y    k    )  ,and use                     y        ^                    k      +      1       as your initial guess for       y          k      +      1      .

Answered question

Answer & Explanation

New Questions in Integral Calculus