I need to solve following ODE

dU/dt=S(U) in which S(U) is geometric source vector in symmetric Euler equations.

Due to instabilities resulting from stiffness, I am considering the use of an implicit time differential scheme. Among them, I first think about the Trapezoidal rule as

U^n+1=U^n+(S(U^n)+S(U^n+1))/dt.

The problem is that I do not know how to construct S(U^n+1). Could you anyone give me idea or direction to find books?

Thanks.

You can expand S(U^n+1) using taylor series. Look up in any CFD book for implicit time marching.

Do you mean that S(U^n+1)=S(U^n)+dS/dU*(U^n+1-U^n)? The problem in my implementations is that S consists of geometric term and any source term to impose sponge layer effects at boundaries. Thus, it was not easy to expand dS/dU for complex relations. Anyway thanks for your reply.

to solve u' = F(u,t), where functionality of F is not explicit in relation with u and t

there areseveral methods that only need time step size and a function that evaluate F at specific "u" and "t", e.g. Rungee Kuta methods (http://en.wikipedia.org/wiki/Runge-Kutta_methods)

if you can not evaluate F at mid way of time step (RK needs), you could take its as picewise constant or use one-sided interpolation from avaialble data.

Did you mention it for implicit schemes? I needed to implement an implicit scheme to continue the computation beyond a certain instant. I am currently suffering from stability problems when the compressible Euler code with a geometric source applied to simulate a wide range of fluid flows.

Thanks in advance!

As you mentioned, I obtained interpolated values in implicit method as

U^n+1=U^n+1/2*dt*(S(u^n)+S(u*^n+1))

where u*^n+1 is interpolated from previous solutions. The procedure to interpolate is based on

U*^n+1=U^n+(U^n-U^n-1). This is my idea to obtain values. Did you intend to let me know this way? Or is there any other way? I really appreciate your great advice

suggested methods are explicit, but with good stability bound, note that the R-K method are like predictor-corrector methods and they try to reconstruct function with higher order and lower needed data, i should say that your extrapolation is also somhow explicit and it is like Adams-Bashforth method which is bit unstable than 4th-R-K, so try R-K methods first, they may resolve. They are extremly simple to implement.

i suggest extrapolation to evaluate f(u,t) at midway. i.e.. at (u+du, t+dt), where u+du<n^{n+1} t+dt<t^{n+1}

what is structure of "F", i do not understand what you mean from "geometric"?

do you have IMSL library, if, it has several prepared routine inter/extrapolation and ODE solvers.

If you know S(u) and if you know U initial value why you dont'use R-K?

I used RK scheme but it caused instabilities when the simulation went long. At early time(within sub-microsecond), it worked well. But at relatively late time(around sub-millisecond), the computation failed to continue after a certain instant. I guessed that this instability results from stiffness of ODE which I was solving. The ODE is a part of compressible Euler equation with a geometric source implementing symmetry spherically or cylindrically. Unfortunately, it leaded to the sudden abortion of computation when I used RK scheme to numerically solve that ODE. It is

dU/dt=S(u) where U=[rho, rho*u, rho*E]^T and S(u)=1/r[rho*u,rho*u^2m rho*(E+p)]^T. This is a simplified version of compressible Euler equation implementing symmetric flow motion. I did not know about this difficulty before. However, when I attempted to see late-time flow motion, I had to run the compressible code somewhat long and I met with sudden abortion of computation due to several reasons. I analyzed that the stiffness of ODE is a main source of this instability. To overcome this issue, I tried to implement an implicit time integration which is free from CFL restriction.

If you have any other advices, could you give me that? I am looking forward to having change to talk all together.

Thanks.

if you need a really implicit solver, you could perform an outer iteration within each time step and go uot when your predicted u^n+1 converge.

there are several implicit ode solver but cost is high, e.g. see Adamsâ€"Moulton methods

http://en.wikipedia.org/wiki/Linear_multistep_method

if you look at a classic numerical analysis book found formal information needed to solve ode

e.g. this one is good: Numerical Analysis, R. Burden

you may find this paper interesting:

http://www.global-sci.com/freedownload/v3_479.pdf