# Solving first order linear ODE numerically

Given a first order linear ODE of the form $$\dot{y} + p(t) y = q(t)$$ with initial condition $$y(\tau_1) = y_0$$, I want to find the value of $$y(\tau_2)$$. What is the most efficient way to do this?

One potential solution that I thought of was to use the solution of linear first order ODE with integrating factor $$y(t) = e^{-h} \left( \int e^h q dt + c \right), \; h = \int p(t) dt$$

and substitute the values of $$\tau_1$$ to find $$C$$ and $$\tau_2$$ into it to find $$y(\tau_2)$$.

But, in my particular problem, the indefinite integral $$h = \int p(t) dt$$ cannot be calculated analytically. And I'm not sure how to handle the $$h$$ numerically. Thoughts?

• Quadrature rules for numerical integration exists (trapezium, Simpson, 3/8). Runge-Kutta methods and exponential Runge-Kutta methods for numerical ODE solutions also exist. ... What is your aim here, a specialized method for first-order linear DE, an exploration of general methods on the example of first-order linear DE, something even more general or specialized? Sep 16, 2022 at 9:39
• I'm particularly interested in specialized methods for first-order linear ODEs. Sep 16, 2022 at 9:43
• The $h$ as it is not well defined. It is unique only up to a constant. You could define $h(t)=\int_0^t p(s) ds$ (or from some other point) and then at each point $t$ evaluate it numerically. You would also have to evaluate the integral in the solution numerically (after plugging in some end values). However evaluating two integrals numerically may be more expensive then a numerical solution to the ODE (although its not that bad). You might be better solving the ODE with some kind of Runge-Kutta method (e.g Euler or Heun).
– user765629
Sep 16, 2022 at 12:24
• Sometimes you can also get away with the power series decompositions. Sep 24, 2022 at 0:50
• What is $\tau_2$ ? Have you any assumption on the periodicity of $p$ of $y$ ? Have you the complete problem ? Is is supposed to be solved numerically as your professor or your book presents to you ?
– EDX
Sep 24, 2022 at 19:56

I agree with many of the comments that the usual way to approach this numerically is just to consider it as the first order equation $$\dot y = f(t,y)$$ (where here $$f(t,y) = q(t) - p(t)y$$) and then use any standard ODE method like Euler or Runge-Kutta. But I feel like you're looking for a method specifically utilizing the linear structure of the problem, so some thoughts:

Even though you cannot integrate $$p(t)$$ analytically, you can still numerically integrate it, for example by using trapezoidal rule or Simpson's rule. Going this route, if you chose $$h(t) := \int_{\tau_1}^tp(s)ds$$, then we have that $$\frac{d}{dt}(e^{h(t)}y(t)) = e^{h(t)}q(t),$$ so solving the problem amounts to numerically computing two integrals, like so:

We are given $$y_0 := y(\tau_1)$$ as an initial value and $$h_0 := h(\tau_1) = 0$$ by construction. Now assume, taking steps of size $$k$$, that we know $$y_i, h_i$$ at time step $$t_i = \tau_1 + ik$$. Then if we use say the trapezoid rule, we have that $$h_{i+1} = h_i + \frac{k}{2}(p(t_i) + p(t_{i+1}))$$ and $$e^{h_{i+1}}y_{i+1} = e^{h_i}y_i + \frac{k}{2}(e^{h_{i}}q(t_i) + e^{h_{i+1}}q(t_{i+1}))$$, which after eliminating the $$h_i$$ parameters turns into the following formula for stepping $$y_i$$ in time directly: $$y_{i+1} = e^{-\Delta}y_i + \frac{k}{2}\left(e^{-\Delta}q(t_i) + q(t_{i+1})\right), \qquad\qquad \Delta = \frac{k}{2}(p(t_i) + p(t_{i+1}))$$ So numerically this formula only involves evaluating $$p$$ and $$q$$ at various time steps, as we've removed any explicit dependence on $$h$$ (of course it's still there implicitly in the $$\Delta$$ term). Then you step in time $$N = \frac{\tau_2-\tau_1}{k}$$ steps until you've reached $$y_N = y(\tau_2)$$.

While in general this "works" in a loose sense, you do ask about an "efficient" method and I cannot speak to how well this would converge or how fast it would run comparatively (I'll that for you to figure out or ask in another question). But I do want to offer an alternative route, which is again to consider the setup $$\dot y = f(t,y)$$ but realize that the linear structure of the problem means that implicit methods like Backwards Euler or Crank Nicholson are actually fairly efficient to implement; for example, Backwards Euler would be $$y_{i+1} = y_i + kf(t_{i+1},y_{i+1}) = y_i + k(q(t_{i+1}) - p(t_{i+1})y_{i+1})$$ $$\implies y_{i+1} = \frac{y_i + kq(t_{i+1})}{1 + kp(t_{i+1})}$$ These methods have proven convergence rates and work very well for stiff problems, much better than explicit methods with larger time steps. Since I do not know the specific equations you are working with I can't say for sure how helpful this would be but something to consider.

Note: If you are looking for a method that solves this analytically and not numerically with time steps, you will also need $$p(t)$$ to be analytically integrated since the unique solution involves this integral. So based on the information given you have to use a time-step method

Let me try exposing how you can use Newton's method, to find the approximate value of $$y$$ at $$\tau_2$$

Knowing the value of $$\tau_1$$ you want to subdivide the interval $$[\tau_1,\tau_2]$$ into small 'acceptable' time 'steps', call that time step $$\Delta t= \frac{\tau_2- \tau_1}{N}$$ and call the instants of time $$t_1, t_2,... , t_k,...$$ , here $$t_1= \tau_1+\Delta t$$ and $$t_k= \tau_1+k\Delta t$$

You start by estimating the value $$y(t_1)$$ by just considering that this value is $$y(\tau_1)\ +$$ whatever a linear progression is allowing, that is $$y(t_1)= y(\tau_1)+ y'(\tau_1)\Delta t= y(\tau_1)+ \Big(q(\tau_1)-p(\tau_1)y(\tau_1)\Big)\Delta t$$ Notice how you use the differential equation to compute the slope.

Continue like this to compute each $$y(t_k)= y(t_{k-1})+ \Big(q(t_{k-1})-p(t_{k-1})y(t_{k-1})\Big)\Delta t$$

You will get more accurate result whenever you take a bigger $$N$$, that is, when you refine your subdivision.

A good way to calibrate would be to notice how much the functions $$p, q$$ vary, it makes sense to consider a time step where the variations of $$p$$ and $$q$$ are both small.