Issues with coupled 1D-ODE

General Discussion
1

Hi there,
I’m rather new with FreeFEM and I’m trying to solve a simple 1D-ODE (2nd-order) with boundary conditions on the same (leftmost) boundary.
Since it is not possible in the weak formulation to impose Dirichlet- and von Neumann-conditions at the same point, I set up a coupled system of 2 ODEs where I can prescribe Dirichlet-data separately for both unknowns. More specifically, let’s consider

u’’ = u - 2sin(x)
u(0) = 0, u’(0) =1,

the obvious solution being u(x) = sin(x).
The code I came up with, is shown here:

// grid definition: we want n points in the interval [0,L]
int n = 1000;
real L = 15*pi;


// specify label indicators for left- and rightmost grid point
int [int] labs = [3, 4];

// create 1D meshL-object
meshL Th = segment(n, [x*L], label = labs);


// create FE-space, test- and trial-functions
fespace Vh(Th, P1);
Vh u1, u2, v1, v2;

// function definition for righthand-side
func f = -2.0*sin(x);


//---------------------------------------------------------------------------
// we want to solve y'' = y - 2*sin(x) = y + f with BCs y(0) = 0; y'(0) = 1
// Dirichlet- and von-Neumann-BC at the same point is not possible in weak 
// formulation, therefore recast into first-order system via
// y1' = y2
// y2' = y1 - 2*sin(x) 
// with pure Dirichlet-BCs y1(0) = 0; y2(0) = 1
//---------------------------------------------------------------------------

// set up and solve weak form
solve simple([u1, u2], [v1, v2], solver = UMFPACK)
    = int1d(Th)(
        dx(u1) * v1
      - u2 * v1 
      + dx(u2) * v2
      - u1 * v2

    )
    - int1d(Th)( 
        f*v2
    )
    + on(3, u1=0.0, u2 = 1.0);

However, when plotting the solution while the left part looks fine and all boundary conditions are fulfilled, the solution gets progressively worse on the right side of the interval and I have no idea why and what I can do to prevent it..

Any help is greatly appreciated, thanks!

Ps: I tried to dig into various example codes and forum posts but wasn’t able to find something that could help me, but I might have overlooked something..

2

Hello,
I tried you code. From my point of view it works rather well.
I plotted the result for n=100

sin
sin681×435 3.5 KB

The black line is the exact solution, the red line is the computed one.
There is only a discrepancy close to the right boundary.

3

Hi François,

thanks for your reply. I agree, at first sight it doesn’t look so bad; although I would have anticipated for the solution to behave perfectly well also at the right boundary. It’s not clear to me why the accuracy goes down there.
I gets worse when you plot u1 and u2 together; u2 should be a pure cosine function and it’s even worse than the behaviour for u1 (see below):

Bildschirmfoto 2025-04-21 um 13.47.01
Bildschirmfoto 2025-04-21 um 13.47.011408×852 69.5 KB

To get the code closer to the expected result, I made some change to the FE-space via

fespace Vh1(Th, P2);
fespace Vh2(Th, P1);
Vh1 u1, v1;
Vh2 u2, v2;

Perhaps surprisingly that leads to perfect agreement between theory and simulation for interval lenghts that are even multiples of 2 pi, while it fails at the right boundary for interval lenghts that are odd multiples of 2 pi.
I’m rather clueless why this happens, in particular since a simple test implementation of the ODE

y’ = cos(x) , y(0) = 0

in FreeFEM has no problem whatsoever to get the correct solution at the right boundary for any intervall lengths..

4

You can use + int1d(Th, qfe=qf1pElump)((abs(x) < 1e-5)*v1); to add a Neumann boundary condition on the left side.

5

in the equation y’=cos(x) there is no growth (the homogeneous equation y’=0 has no solution with exponential growth). Moreover this equation is only first-order.
For y’'=y-2sin(x), the homogeneous equation has exp(x) as solution. This implies difficulty to have a stable scheme.
In order to get a “good” scheme we would need:
-stability
-accuracy at the right boundary (the discretization should be consistent with the equation at sufficiently high order), which is probably missing here

6

Hi Quentin,
thanks for your suggestion. I tried to modify the code and used this snippet

// set up and solve weak form
solve simple(u, v, solver = UMFPACK)
    = int1d(Th)(
      - dx(u) * v
      - u * v 
    )
    - int1d(Th)( 
        f*v
    )
    + int1d(Th, qfe=qf1pElump)(
        (abs(x) < 1e-5)*v
    )
    + on(3, u=0.0)
    ;

as a replacement but it didn’t seem to properly compute the analytic solution (see below):

Bildschirmfoto 2025-04-21 um 17.09.52
Bildschirmfoto 2025-04-21 um 17.09.521408×866 67.6 KB

7

maybe you need to change the variational form, take the integral by part to u’'v, to have the correct boundary term.

8

Certainly this has to to with stability at one point or the other. It’s still strange to me that the code leads to the correct solution for “properly chosen” intervals. But maybe this is due to the reason that the correct solution has a root there and this helps for convergence.
I’m a little bit surprised that the FEM-formulation is not as trivial as I expected, especially when compared to simple finite difference-approaches but maybe FEM is not the best tool for IVPs. I initially thought, my little toy problem would be rather easy to implement…

9

I think it is, not sure

solve simple(u, v, solver = UMFPACK)
    = int1d(Th)(
      - dx(u) * dx(v)
      - u * v 
    )
    - int1d(Th)( 
        f*v
    )
    + int1d(Th, qfe=qf1pElump)(
        (abs(x) < 1e-5)*v
    )
    + on(3, u=0.0)
    ;

but here we miss the right Neumann condition, it perhaps will give u’=0 on the right boundary.

10

Jep, this is indeed what is happening :grinning_face:

11

The fact is that your problem is a Cauchy problem with all boundary conditions on the same side. It is not a boundary value problem. This is why a variational formulation is not the best way to describe it. The usual treatment is via an ODE solver for which known (u1_n,u2_n) it computes (u1_{n+1},u2_{n+1}).
The finite element method does not preserve the evolution structure of the continuous problem: the solution at time t (thinking of the position x as a time is more intuitive) should only depend on the past. Here with the finite element method, if you compute u on [0,L], the restriction of u on [0,L/2] is not the solution to the variational formulation on the interval [0,L/2].

12

Thanks for this clarifying remarks; however I still hope to find some literature concerning IVPs and the FE-method, I’m not giving up this fast :laughing: