# implicit odes solution using fdm

I am solving a non linear second order implicit initial value problem using finite difference method, but my results do not converge. Please guide me with an example, how we can apply finite difference method with big accuracy. The problem is $$(A+By(t)-C(-y'(t))^-1/3)y''(t)=-(By'(t)^2+C),\qquad y(t=0)=0,$$ using fdm, time can be from zero to 2 mints. A,B and C are constants.

• I guess, this question cannot be answered in general. You should provide more information about the type of problem you are trying to solve. – H. Rittich Aug 4 '17 at 13:12
• If you have access to matlab, you can solve this with the ode15i function, mathworks.com/help/matlab/ref/ode15i.html. If you don't have matlab, the method used in ode15i is discussed in this paper by Shampine, faculty.smu.edu/shampine/cic.pdf. – Bill Greene Aug 4 '17 at 14:13
• Welcome to Scicomp.SE! You should add more details before we can help you: Which finite difference method are you using (there are uncountably many...)? What exactly do you mean by "not converging"? – Christian Clason Aug 4 '17 at 14:47
• I want to apply backward difference method (implicit) in time and then proceeding through iterations in matlab. Converging means me that, i applied explicit method and i could not succeeded to achieve results in paper (research article). – B. Munir Aug 5 '17 at 5:04

You also have to specify which scheme you want to use.

Let's say that you want to use the simple Backward-Euler (BE) scheme. For BE scheme, your governing equation becomes

$\left(A+By_{n+1}-C(y'_{n+1})^{-1/3}\right)y''_{n+1}=-\left(B(y'_{n+1})^2+C\right) \tag{1}$ with $y'_{n+1}=\frac{y_{n+1}-y_n}{\Delta t} \tag{2}$

$y''_{n+1}=\frac{y'_{n+1}-y'_n}{\Delta t} \tag{3}$

Using (2) in (3), $y''_{n+1}=\frac{1}{\Delta t^2}[y_{n+1}-y_n] - \frac{y'_n}{\Delta t} \tag{4}$

Now, using (2) and (4) in (1), you will get a nonlinear equation of the form

$f(y_n,y'_n, y_{n+1})=0 \tag{5}$ which can be solved for $y_{n+1}$ using an iterative scheme, for example, Newton-Raphson scheme. Once you get $y_{n+1}$ you can compute velocity and acceleration from (2) and (4), respectively.

You can follow the same procedure for other schemes. Equations (1), (2), (3) and (4) will be different for different schemes.

Since you have a second order ODE you must provide 2 initial conditions, say: $$y'(0) = y'_0 \qquad y(0) = y_0\tag{1}$$

One thing you can do is transform the original second order equation: $$\left(A+By-Cy'^{-1/3}\right)y''=-\left(By'^2+C\right) \tag{2}$$ into a system of first order ODEs $$\left\{\begin{array}{ll}% v' =& -\frac{Bv^2+C}{A+By-Cv^{-1/3}}\\ y'=&v \end{array}\right. \tag{3}$$ with initial conditions: $$v(0) = y'_0 \qquad y(0) = y_0$$

System $(2)$ can be put into one general vector equation: being $\vec{x} = [x_1,x_2]^{T}$ for $x_1=v$ and $x_2=y$ and $\vec{f} = \left[-\frac{Bv^2+C}{A+By-Cv^{-1/3}},\,v\right]^T$. Therefore: $$\frac{d \vec{x}}{dt} = \vec{f}(\vec{x})$$

Yu can solve it easily with matlab using any ODE solver that supports stiff equations (implicit numerical method).

• thank you but I am interested in FDM, as said above again. – B. Munir Aug 5 '17 at 5:05