The answer to your question requires a little of explanation.

There are several different Finite Element Methods, not just one. ALL of these methods can be classified in 2 groups: (1) F.E. methods based on fundamental Lemma of Calculus of Variations (Galerkin method, Petrov-Galerkin, Weighted residual, Galerkin method with weak form) and (2) F.E. method based on a residual functional associated with the differential equation (least-squares method).

There are several types of differential equations. Hyberbolic, Elliptic, Parabolic, First order, Second order, etc. HOWEVER, ALL differential equations have a differential operator associated with them and ALL differential operators can be classified in 3 groups: (1) self-adjoint (2) non-self adjoint (3) non-linear.

For a F.E. method to be used successfully for a specific differential equation, the resulting coefficient matrix in [K]{u}={F} must be positive-definite (This means that [K] can be inverted and the obtained solution {u} will be unique).

The following 3 statements are not opinions but Theorems with proofs:

- When the differential operator is self-adjoint: Galerkin method with weak form and least squares process are the only F.E. methods that guarantee positive-definite coefficient matrix [K].

- When the differential operator is non-self adjoint or non-linear: least square process is the only F.E. method that guarantees a positive-definite coefficient matrix [K].

- When considering an IVP, the differential operator is either non-self adjoint or nonlinear. It is never self-adjoint. Hence, the only F.E. method that can guarantee positive definite coefficient matrices for all IVPs is least squares process.

There is a lot of research done to make the Galerkin method with weak form produce coefficient matrices that are positive-definite. These are known as stabilizing techniques. However, all of these techniques, without exception, end up changing the original IVP.

Many people are under the impression that Galerkin method with weak form is the only F.E. method. This is because ALL commercial F.E. softwares are based on this method. The reason for this being that the differential operators of differential equations describing elastic solid mechanics are always self-adjoint, hence Galerkin method with weak form works great! Least-squares process also works great, however, Least-squares process requires additional resources (in terms of interpolation theory) that didn't exist at the time when commercial F.E. softwares started to appear (1960s). Hence, people didn't consider Least-squares process. Around the late 1990s, those additional resources needed for least-squares were invented, and least-squares started being used in the research area. However, since so many companies were built on Galerkin  method with weak form and they are only marketed towards solid mechanics applications, they still continue to use it.

Other methods to solve differential equations such as: Finite Volume, Finite Difference, Runge-Kutta etc suffer from the same disease than Galerkin method with weak form. Meaning, that they do not guarantee positive-definite coefficient matrices. This does not mean that unique solutions are not possible using these methods, but it means that you will only know it after the [K]{u}={F} is obtained and if [K] is not positive definite, then stabilizing techniques must be used. Runge-Kutta method and Finite Difference methods are very easy and fast to implement. Hence they popularity. However,  when considering complex problems where you might end up having a [K] matrix that is a million-by-million, then they might or might not work and if they don't work, stabilizing techniques might or might not give you the right solution.

So to answer your question, F.E. method can be used to solve IVP. Either using Galerkin method with weak form combined with stabilizing techniques (see works of T. Hughes) or using least-squares process. If you are interested on using F.E. method to solve ODEs in time, then you should real this article:

- K. S. Surana, L. Euler, J. N. Reddy, A. Romkes: Methods of Approximation in hpk Framework for ODEs in Time Resulting from Decoupling of Space and Time in IVPs. American J. Computational Mathematics 1(2): 83-103 (2011)

This article explains into details how to handle IVPs by decoupling space and time and as a result, obtaining a system of ODEs in time. Then it compares F.E. methods and other methods such as Newmark, Wilson method. I hope that gives you a better idea of the subject.

For initial value problems conditions are given at a single point t=0 and t >0, so Taylor's expansion is possible and a solution can be obtained in the near t=0. Continuing this way solutions can be obtained. However, for BVP the conditions are given at two different points so either a shooting method or FDM or FEM is to be applied.

   For IVP , IVP can not be applied because an integral formulation can not be obtained to apply FEM.

The answer to your question requires a little of explanation.

There are several different Finite Element Methods, not just one. ALL of these methods can be classified in 2 groups: (1) F.E. methods based on fundamental Lemma of Calculus of Variations (Galerkin method, Petrov-Galerkin, Weighted residual, Galerkin method with weak form) and (2) F.E. method based on a residual functional associated with the differential equation (least-squares method).

There are several types of differential equations. Hyberbolic, Elliptic, Parabolic, First order, Second order, etc. HOWEVER, ALL differential equations have a differential operator associated with them and ALL differential operators can be classified in 3 groups: (1) self-adjoint (2) non-self adjoint (3) non-linear.

For a F.E. method to be used successfully for a specific differential equation, the resulting coefficient matrix in [K]{u}={F} must be positive-definite (This means that [K] can be inverted and the obtained solution {u} will be unique).

The following 3 statements are not opinions but Theorems with proofs:

- When the differential operator is self-adjoint: Galerkin method with weak form and least squares process are the only F.E. methods that guarantee positive-definite coefficient matrix [K].

- When the differential operator is non-self adjoint or non-linear: least square process is the only F.E. method that guarantees a positive-definite coefficient matrix [K].

- When considering an IVP, the differential operator is either non-self adjoint or nonlinear. It is never self-adjoint. Hence, the only F.E. method that can guarantee positive definite coefficient matrices for all IVPs is least squares process.

There is a lot of research done to make the Galerkin method with weak form produce coefficient matrices that are positive-definite. These are known as stabilizing techniques. However, all of these techniques, without exception, end up changing the original IVP.

Many people are under the impression that Galerkin method with weak form is the only F.E. method. This is because ALL commercial F.E. softwares are based on this method. The reason for this being that the differential operators of differential equations describing elastic solid mechanics are always self-adjoint, hence Galerkin method with weak form works great! Least-squares process also works great, however, Least-squares process requires additional resources (in terms of interpolation theory) that didn't exist at the time when commercial F.E. softwares started to appear (1960s). Hence, people didn't consider Least-squares process. Around the late 1990s, those additional resources needed for least-squares were invented, and least-squares started being used in the research area. However, since so many companies were built on Galerkin  method with weak form and they are only marketed towards solid mechanics applications, they still continue to use it.

Other methods to solve differential equations such as: Finite Volume, Finite Difference, Runge-Kutta etc suffer from the same disease than Galerkin method with weak form. Meaning, that they do not guarantee positive-definite coefficient matrices. This does not mean that unique solutions are not possible using these methods, but it means that you will only know it after the [K]{u}={F} is obtained and if [K] is not positive definite, then stabilizing techniques must be used. Runge-Kutta method and Finite Difference methods are very easy and fast to implement. Hence they popularity. However,  when considering complex problems where you might end up having a [K] matrix that is a million-by-million, then they might or might not work and if they don't work, stabilizing techniques might or might not give you the right solution.

So to answer your question, F.E. method can be used to solve IVP. Either using Galerkin method with weak form combined with stabilizing techniques (see works of T. Hughes) or using least-squares process. If you are interested on using F.E. method to solve ODEs in time, then you should real this article:

- K. S. Surana, L. Euler, J. N. Reddy, A. Romkes: Methods of Approximation in hpk Framework for ODEs in Time Resulting from Decoupling of Space and Time in IVPs. American J. Computational Mathematics 1(2): 83-103 (2011)

This article explains into details how to handle IVPs by decoupling space and time and as a result, obtaining a system of ODEs in time. Then it compares F.E. methods and other methods such as Newmark, Wilson method. I hope that gives you a better idea of the subject.

@Daniel : I do not agree with you on the following statement

When considering an IVP, the differential operator is either non-self adjoint or nonlinear. It is never self-adjoint. Hence, the only F.E. method that can guarantee positive definite coefficient matrices for all IVPs is least squares process.

  Whether a differential equation is a BVP or IVP depends on the conditions for existence of its unique solutions given. It is nothing to do with differential operator. Same differential equation can be a BVP or IVP depending on if the conditions are given at one point or at two points ( for second order ODE). The equation may be in self adjoint or linear form yet it is IVP. IVP always have stability issue .

The basic points of a FEM:

1. Define the trial function space

2. define a suitable condition to determine the coefficient of the trial functions expansion that is the best approximation of the solution in the trial functional space.

From this point of view the answer to your question is :

NO the FEM can be used to approximate in space and time the unknown function but it is not always effective respect to other approximating methods.

It is possible to use FEM for IVPs (ODEs), but methods based on FD are popular because of the simplicity in implemention. The approximate solution using FEM is stable and high-order accurate provided that the continuous problem itself is stable. If you consider FEM for such problems here is a brief outline of how it is implemented:

   (The continuous problem)                     Du = 0;     with u(0) = u_0

1. Obviously, discretise the time domain (i.e. the domain of the problem say [0, L]).

    0 = t_0 < t_1 < . . . < t_N = L; with h_i =  t_i+1 - t_i. 

2. Start from the first subinterval [x_0, x_1]. Then enforce the problem weakly (Galerkin method) on the interval using finite number of basis functions (in FEM terminology-shape functions)  and use the initial condition as the boundary value (Dirichlet boundary condition) and solve the matrix the matrix equation, now you get the solution at time t_1 and use it as initial condition (Dirichlet boundary condition) to solve the problem on the next subinterval [x_1, x_2] by applying the same procedure. Continue this way until you reach the desired time stage.

Note that there are several variation of FEM for IVPs depending of the necessity of stabilization; to mention a couple of them --time-discontinuous DG FEM and Least-square FEM.

Please, my question is restricted to ODES both linear and nonlinear.Consider the 1st order nonlinear ODE ux + u2 = 0. with the initial value u(0) = 0. The exact solution is u = 1/(1+x). This problem can be solved numerically via Galerkin FEM without any stabilization. My question is can the IVP for the Bessel ODE xu'' + u' + xu = 0 given u(0) = 1.00 and u'(0) = 0 be solved using FEM. It is straight forward to solve it using the Runge-Kutta methods! On the other hand the BVP xu'' + u' + xu = 0 given u(0) = 1.00 and u(1) = 0.7652 can be  be solved readily using the FEM.

Yes, the second order initial value problem can be solved using FEM. You can still write it as a system of first order DE, like

              u' = w

              xw'=w + xu   with u(0) = 1. and w(0) = 0.

  Then it can be solved using the FEM procedure given in my previous comment, .

Please Mebratu

  Using your suggestion we shall be working with two unknowns u and w. Which implies a cubic interpolation function is the minimum. Could you please give me the 'stiffness matrix' for the solution of Bessel eqn using your suggestion. I hope I am not asking too much?

Dear Amaechi,

First I don't agree with the minimum cubic interpolation function. The in the following FE implementation I used linear shape functions for both u and v. The residula --R-- and Tangent (Stiffness)--K-- are generated using Matlab's symbolic tool

R = [(-u1/2 + u2/2 -dx*w1/6 - dx*w2/3);...

        (dx^2*u1/12 + dx^2*u2/4 - dx*w1/6 + 2*dx*w2/3 +dx*u1*xn/6 +...

         dx*u2*xn/3 - w1*xn/2 + w2*xn/2)];

K = [          1/2,                               -dx/3;

         dx^2/4 + dx*xn/3,        2*dx/3 + xn/2];

Where dx is step length i.e dx = x_n - x_n-1;

             u1, w1 are the previous  step values at x_n-1;

             u2, w2 are the unknowns at the current step x_n;

            xn the previous step i.e xn = x_n;

Note the nonlinearity is going to be treated by 2 to 3 Newton's iteration steps.

Hope this can clarify the stuff; See the attached plot of first Bessel's function approximation by this method. Please Note the following

I  have  also encountered the ivp s solved by FEM. For example,  the  following review paper by Cockburn

https://www.ima.umn.edu/preprints/may2003/1921.pdf

solves a first order ivp by discontionus Galerkin FEM .   It seems to me that the proposed approach can be extented to second order ivp's.  

However, many standard textbooks do  not  mention the applicability of FEM methods to ivp's. Although I  am not   a practitioner  of FEM  methods, as far as  I know    FEM is  especially preferred  due to its ability  to handle  the problems with complex geometry in spaces with desired accuracy.  But,  time  geometry  consists of only  one sided  line.  Very  simple to  make partition  (use adaptivity if  necessary)   and   march in time  (implicit or explicit).

As a  summary, yes  FEM  can  be used BUT    should it  be used  in practice ?   

Hello Mebratu Wakeni

Thanks for your valuable answers. I don't understand how the R and K matrices are used to obtain the solution of Bessel equation. My computer cannot open the".eps" you provided. Could you supply another file format like matlab ".m" or a text file.

Hi Amaechi

Please find attached the Matlab script file.

Best Regards,

Mebratu

Thank you very much my friend

Mebratu Wakeni · I appreciate your cooperation in providing me with the Matlab script file.  

However, I have a question. Why is it that your stiffness matrix is a 2 by 2 matrix instead of a 4 by 4 matrix since we are dealing with two unknowns  w and u. If we use linear interpolation the unknowns in each element should be w1, w2, u1, u2 ?

Regards

Amaechi