# Issues with coupled 1D-ODE **URL:** \ **Category:** General Discussion\ **Created:** [April 20, 2025, 8:50am UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872 "2025-04-20T08:50:30Z")\ **Posts on this page:** 12\ **Page:** 1
**Author:** ![erigae](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/erigae/32/3130_2.png) [@erigae](https://community.freefem.org/u/erigae)\ **Post date:** [April 20, 2025, 8:50am UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/1 "2025-04-20T08:50:30Z")
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: ```auto // 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.. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [April 21, 2025, 11:03am UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/2 "2025-04-21T11:03:02Z")
Hello, I tried you code. From my point of view it works rather well. I plotted the result for n=100 ![sin](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/b/b7b0067bbe49269e69c0058c0c64a850f42991ec.png) The black line is the exact solution, the red line is the computed one. There is only a discrepancy close to the right boundary. ---
**Author:** ![erigae](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/erigae/32/3130_2.png) [@erigae](https://community.freefem.org/u/erigae)\ **Post date:** [April 21, 2025, 12:07pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/3 "2025-04-21T12:07:51Z")
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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/36e7efa5f0506ea26913227741ef258a35ff37bf.png) To get the code closer to the expected result, I made some change to the FE-space via ```auto 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.. ---
**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\ **Post date:** [April 21, 2025, 1:31pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/4 "2025-04-21T13:31:19Z")
You can use` + int1d(Th, qfe=qf1pElump)((abs(x) < 1e-5)*v1);` to add a Neumann boundary condition on the left side. ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [April 21, 2025, 2:28pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/5 "2025-04-21T14:28:03Z")
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 ---
**Author:** ![erigae](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/erigae/32/3130_2.png) [@erigae](https://community.freefem.org/u/erigae)\ **Post date:** [April 21, 2025, 3:10pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/6 "2025-04-21T15:10:24Z")
Hi Quentin, thanks for your suggestion. I tried to modify the code and used this snippet ```auto // 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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/6787dd968e10f4169884dcb557631017d2bd74cc.png) ---
**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\ **Post date:** [April 21, 2025, 3:15pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/7 "2025-04-21T15:15:50Z")
maybe you need to change the variational form, take the integral by part to u’'v, to have the correct boundary term. ---
**Author:** ![erigae](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/erigae/32/3130_2.png) [@erigae](https://community.freefem.org/u/erigae)\ **Post date:** [April 21, 2025, 3:16pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/8 "2025-04-21T15:16:21Z")
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… ---
**Author:** ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)\ **Post date:** [April 21, 2025, 3:19pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/9 "2025-04-21T15:19:26Z")
I think it is, not sure ```auto 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. ---
**Author:** ![erigae](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/erigae/32/3130_2.png) [@erigae](https://community.freefem.org/u/erigae)\ **Post date:** [April 21, 2025, 3:23pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/10 "2025-04-21T15:23:50Z")
Jep, this is indeed what is happening 😀 ---
**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\ **Post date:** [April 22, 2025, 1:41pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/11 "2025-04-22T13:41:40Z")
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]. ---
**Author:** ![erigae](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/erigae/32/3130_2.png) [@erigae](https://community.freefem.org/u/erigae)\ **Post date:** [April 22, 2025, 6:30pm UTC](https://community.freefem.org/t/issues-with-coupled-1d-ode/3872/12 "2025-04-22T18:30:20Z")
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 😆