# Solving the Fractional Differential Equation $\alpha D^{2v}y(x) + \beta D^{v}y(x) + \gamma y(x) = f(x)$

### Introduction

I'm very interested in fractional calculus, especially in Fractional Differential Equations (FDEs). The question arises how to solve the FDE $$\alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + \gamma \cdot y\left( x \right) = f\left( x \right)$$ where $$\left\{ \alpha,\, \beta,\, \gamma,\, y\left( x \right),\, f\left( x \right),\, v,\, x \right\} \in \mathbb{C}$$ and where $$\operatorname{D}^{v} \left( y\left( x \right) \right)$$ is the $$v$$-th derivative of $$y\left( x \right)$$. The FDE is simple, but I don't see a trivial solution. My attempts at a solution so far have only led to a lesson or restricted the FDE to an ODE with $$v \in \mathbb{R}$$. Hence the question: How to solve the FDE $$\alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + \gamma \cdot y\left( x \right) = f\left( x \right)$$ where $$\left\{ \alpha,\, \beta,\, \gamma,\, y\left( x \right),\, f\left( x \right),\, v,\, x \right\} \in \mathbb{C}$$?

### Similar FDE

Of course there are similar FDEs, like the special case $$f\left( x \right) = 0$$: \begin{align*} \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + \gamma \cdot y\left( x \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + \gamma \cdot y\left( x \right) &= 0\\ \end{align*}

\begin{align*} \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + \gamma \cdot y\left( x \right) &= f\left( x \right)\\ \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \frac{\beta}{\alpha} \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + \frac{\gamma}{\alpha} \cdot y\left( x \right) &= \frac{1}{\alpha} \cdot f\left( x \right)\\ \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + a \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + b \cdot y\left( x \right) &= \frac{1}{\alpha} \cdot f\left( x \right)\\ \end{align*}

\begin{align*} \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \frac{\beta}{\alpha} \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + \frac{\gamma}{\alpha} \cdot y\left( x \right) &= 0\\ \end{align*}

This FDE has e.g. given a solution by: \begin{align*} y\left( x \right) &= \begin{cases} \sum\limits_{k = 1}^{\frac{1}{v} - 1}\left( a^{\frac{1}{v} - k - 1} \cdot \operatorname{E}_{x}\left( -k \cdot v;\, a^{\frac{1}{v}} \right) - b^{\frac{1}{v} - k - 1} \cdot \operatorname{E}_{x}\left( -k \cdot v;\, b^{\frac{1}{v}} \right) \right), &\text{if } a \ne b\\ x \cdot \exp\left( a \cdot x \right) \cdot \sum\limits_{k = -\frac{1}{v} + 1}^{\frac{1}{v} - 1}\left( a^{k} \cdot \left( \frac{1}{v} - \left| k \right| \right) \cdot \exp\left( a^{\frac{1}{v}} \cdot x \right) \right), &\text{if } a = b \ne 0\\ \frac{1}{\Gamma\left( 2 \cdot v \right)} \cdot x^{2 \cdot v - 1 }, &\text{if } a = b = 0\\ \end{cases}\\ \end{align*}

where $$\Gamma\left( \cdot \right)$$ is the Gamma Function and $$\operatorname{E}_{t}\left( \cdot;\, \cdot \right)$$ is the $$\operatorname{E}_{t}$$ Function.

### My Best Attempts At Finding A Solution

#### Frobenius Method

##### Without Substitution

I would assume that $$y\left( x \right)$$ has a Maclaurin Series given by $$y\left( x \right) = \sum\limits_{k = 0}^{\infty}\left( y_{k} \cdot x^{k} \right)$$. Then we could use Euler's approach for the Fractional Derivative of power functions $$\operatorname{D}^{n}\left( x^{m} \right) = \frac{\Gamma\left( m + 1 \right)}{\Gamma\left( m - n + 1 \right)} \cdot x^{m - n}$$ and get this: \begin{align*} \operatorname{D}^{n} \left( y\left( x \right) \right) &= \operatorname{D}^{n} \left( \sum\limits_{k = 0}^{\infty}\left( y_{k} \cdot x^{k} \right) \right)\\ \operatorname{D}^{n} \left( y\left( x \right) \right) &= \sum\limits_{k = 0}^{\infty}\left( y_{k} \cdot \frac{\Gamma\left( k + 1 \right)}{\Gamma\left( k - n + 1 \right)} \cdot x^{k - n} \right)\\ \end{align*}

If we plug $$\operatorname{D}^{n} \left( y\left( x \right) \right) = \sum\limits_{k = 0}^{\infty}\left( y_{k} \cdot \frac{\Gamma\left( k + 1 \right)}{\Gamma\left( k - n + 1 \right)} \cdot x^{k - n} \right)$$ in, we get: \begin{align*} \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + a \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + b \cdot y\left( x \right) &= \frac{1}{\alpha} \cdot f\left( x \right)\\ \sum\limits_{k = 0}^{\infty}\left( y_{k} \cdot \Gamma\left( k + 1 \right) \cdot \left( \frac{1}{\Gamma\left( k - 2 \cdot v + 1 \right)} \cdot x^{k - 2 \cdot v} + a \cdot \frac{1}{\Gamma\left( k - v + 1 \right)} \cdot x^{k - v} + \frac{1}{\Gamma\left( k + 1 \right)} \cdot x^{k} \right) \right) &= \frac{1}{\alpha} \cdot f\left( x \right)\\ \sum\limits_{k = 0}^{\infty}\left( y_{k} \cdot \Gamma\left( k + 1 \right) \cdot \left( \frac{1}{\left( k - 2 \cdot v - 1 \right) \cdot \Gamma\left( k - 2 \cdot v \right)} \cdot x^{k - 2 \cdot v} + a \cdot \frac{1}{\left( k - v - 1 \right) \cdot \Gamma\left( k - v \right)} \cdot x^{k - v} + \frac{1}{\Gamma\left( k + 1 \right)} \cdot x^{k} \right) \right) &= \frac{1}{\alpha} \cdot f\left( x \right)\\ \end{align*}

Where I have no idea how to continue in a meaningful way. Of course, I could substitute the exponents of $$x$$ the same, then do a index-shift to put it back together, hoping that helps, but I'd rather give up on that. However, that would only allow one solution for $$v \in \mathbb{Z}$$, but then it is no longer an FDE, but an ODE.

##### With Substitution (for $$\gamma = 0$$)

The first sensible substitution that comes to mind is $$z\left( x \right) ~{:=}~ \operatorname{D}^{v} \left( y\left( x \right) \right)$$. Which we can rearrange with something and apply the condition $$\operatorname{D}^{a} \left( \operatorname{D}^{b} \left( y\left( x \right) \right) \right) = \operatorname{D}^{a + b} \left( y\left( x \right) \right)$$ so we get this: \begin{align*} \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + 0 \cdot y\left( x \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + 0 &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{v + v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{v} \left( \operatorname{D}^{v} \left( y\left( x \right) \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{v} \left( z\left( x \right) \right) + \beta \cdot z\left( x \right) &= f\left( x \right)\\ \end{align*}

I would assume that $$z\left( x \right)$$ has a Maclaurin Series given by $$z\left( x \right) = \sum\limits_{k = 0}^{\infty}\left( z_{k} \cdot x^{k} \right)$$. Then we could use Euler's approach for the Fractional Derivative of power functions $$\operatorname{D}^{n}\left( x^{m} \right) = \frac{\Gamma\left( m + 1 \right)}{\Gamma\left( m - n + 1 \right)} \cdot x^{m - n}$$ and get this: \begin{align*} \alpha \cdot \operatorname{D}^{v} \left( z\left( x \right) \right) + \beta \cdot z\left( x \right) &= f\left( x \right)\\ \alpha \cdot \sum\limits_{k = 0}^{\infty}\left( z_{k} \cdot \frac{\Gamma\left( k + 1 \right)}{\Gamma\left( k - v + 1 \right)} \cdot x^{k - v} \right) + \beta \cdot \sum\limits_{k = 0}^{\infty}\left( z_{k} \cdot x^{k} \right) &= f\left( x \right)\\ \alpha \cdot \sum\limits_{k = -v}^{\infty}\left( z_{k + v} \cdot \frac{\Gamma\left( k + 1 \right)}{\Gamma\left( k + 1 \right)} \cdot x^{k} \right) + \beta \cdot \sum\limits_{k = 0}^{\infty}\left( z_{k} \cdot x^{k} \right) &= f\left( x \right)\\ \alpha \cdot \sum\limits_{k = -v}^{0}\left( z_{k + v} \cdot \frac{\Gamma\left( k + 1 \right)}{\Gamma\left( k + 1 \right)} \cdot x^{k} \right) + \alpha \cdot \sum\limits_{k = 0}^{\infty}\left( z_{k + v} \cdot \frac{\Gamma\left( k + 1 \right)}{\Gamma\left( k + 1 \right)} \cdot x^{k} \right) + \beta \cdot \sum\limits_{k = 0}^{\infty}\left( z_{k} \cdot x^{k} \right) &= f\left( x \right)\\ \alpha \cdot \sum\limits_{k = -v}^{0}\left( z_{k + v} \cdot \frac{\Gamma\left( k + 1 \right)}{\Gamma\left( k + 1 \right)} \cdot x^{k} \right) + \sum\limits_{k = 0}^{\infty}\left( \left( \alpha \cdot z_{k + v} \cdot \frac{\Gamma\left( k + 1 \right)}{\Gamma\left( k + 1 \right)} + \beta \cdot z_{k} \right) \cdot x^{k} \right) &= f\left( x \right)\\ \end{align*}

However, that would only allow one solution for $$v \in \mathbb{Z}$$, but then it is no longer an FDE, but an ODE.

#### Solving For The Homogeneous Equation

##### With Substitution (for $$\gamma = 0$$)

The first sensible substitution that comes to mind is $$z\left( x \right) ~{:=}~ \operatorname{D}^{v} \left( y\left( x \right) \right)$$. Which we can rearrange with something and apply the condition $$\operatorname{D}^{a} \left( \operatorname{D}^{b} \left( y\left( x \right) \right) \right) = \operatorname{D}^{a + b} \left( y\left( x \right) \right)$$ so we get this: \begin{align*} \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + 0 \cdot y\left( x \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + 0 &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{v + v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{v} \left( \operatorname{D}^{v} \left( y\left( x \right) \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) &= f\left( x \right)\\ \alpha \cdot \operatorname{D}^{v} \left( z\left( x \right) \right) + \beta \cdot z\left( x \right) &= f\left( x \right)\\ \end{align*}

Composing the characteristic equation gives: \begin{align*} \alpha \cdot \operatorname{D}^{v} \left( z\left( x \right) \right) + \beta \cdot z\left( x \right) &= f\left( x \right)\\ \alpha \cdot \lambda^{v} + \beta \cdot 1 &= 0\\ \alpha \cdot \lambda^{v} + \beta &= 0\\ \lambda^{v} &= -\frac{\beta}{\alpha}\\ \lambda &= \left( -\frac{\beta}{\alpha} \right)^{\frac{1}{v}}\\ \lambda_{k} &= \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \operatorname{cis}\left( \operatorname{Arg}\left( -\frac{\beta}{\alpha} \right) + 2 \cdot k \cdot \pi \right)^{\frac{1}{v}}\\ \lambda_{k} &= \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \operatorname{cis}\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right) + 2 \cdot k \cdot \pi}{v} \right)\\ \lambda_{k} &= \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \operatorname{cis}\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right)}{v} + \frac{2}{v} \cdot k \cdot \pi \right)\\ \end{align*}

where $$\operatorname{Arg}\left( z \right)$$ is the argument of z wich is the closest to $$0$$.

This means that the roots of the equation are $$\lambda_{k} = \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \cos\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right)}{v} + \frac{2}{v} \cdot k \cdot \pi\right) + \sqrt[v]{ \left| \frac{\beta}{\alpha} \right|} \cdot \sin\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right)}{v} + \frac{2}{v} \cdot k \cdot \pi\right) \cdot i$$ aka $$\lambda^{v} + \frac{\beta}{\alpha} = \prod\limits_{k = 1}^{\left| \mathbb{S} \right|}\left( \lambda - \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \operatorname{cis}\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right)}{v} + \frac{2}{v} \cdot k \cdot \pi\right) \right)$$ where $$\mathbb{S}$$ is the set of all roots of the equation. This gives the multiplicity of the all roots $$= 1$$ for $$\alpha \ne 0 \ne \beta$$.

With $$z_{h}\left( x \right) = Q_{0}\left( x \right) \cdot \exp\left( \Re\left( \lambda_{k} \right) \cdot x \right) \cdot \cos\left( \Im\left( \lambda_{k} \right) \cdot x \right) + P_{0}\left( x \right) \cdot \exp\left( \Re\left( \lambda_{k} \right) \cdot x \right) \cdot \sin\left( \Im\left( \lambda_{k} \right) \cdot x \right)$$ we'll get \begin{align*} z_{h}\left( x \right) &= c_{2 \cdot k} \cdot \exp\left( \Re\left( \lambda_{k} \right) \cdot x \right) \cdot \cos\left( \Im\left( \lambda_{k} \right) \cdot x \right) + c_{2 \cdot k + 1} \cdot \exp\left( \Re\left( \lambda_{k} \right) \cdot x \right) \cdot \sin\left( \Im\left( \lambda_{k} \right) \cdot x \right)\\ z_{h}\left( x \right) &= \sum\limits_{k = 1}^{\left| S \right|}\left( c_{2 \cdot k} \cdot \exp\left( \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \cos\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right)}{v} + \frac{2}{v} \cdot k \cdot \pi\right) \cdot x \right) \cdot \cos\left( \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \sin\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right)}{v} + \frac{2}{v} \cdot k \cdot \pi\right) \cdot x \right) + c_{2 \cdot k + 1} \cdot \exp\left( \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \cos\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right)}{v} + \frac{2}{v} \cdot k \cdot \pi\right) \cdot x \right) \cdot \sin\left( \sqrt[v]{\left| \frac{\beta}{\alpha} \right|} \cdot \sin\left( \frac{\operatorname{Arg}\left( -\frac{\beta}{\alpha} \right)}{v} + \frac{2}{v} \cdot k \cdot \pi\right) \cdot x \right) \right) \end{align*}

where $$c$$ is a constant. Aka if $$v \in \mathbb{N} \Rightarrow \left| \mathbb{S} \right| = v$$.

Via using $$z\left( x \right) = z_{h}\left( x \right) + z_{p}\left( x \right)$$ we would get: $$\fbox{ y\left( x \right) = \operatorname{D}^{v} \left( z_{h}\left( x \right) + z_{p}\left( x \right) \right) }$$

But I'm not sure how to determine $$z_{p}$$ (Particular Solution), whether I have correctly determined $$z_{h}$$ (Homogeneous Solution) and how to proceed after further determination.

• Are interested in numerical solutions? Apr 5, 2023 at 17:06
• @gciriani I would also be very interested in a numerical solution. If it stays the best answer under my question in the next few days it will get the bounty as well. Apr 5, 2023 at 21:22

You can numerically solve fractional differential equations by using the open-source, Matlab-compatible code that I discussed and used in my answer to a different question: Garrappa, Roberto, Numerical solution of fractional differential equations: a survey and a software tutorial, ZBL06916890.

The software solves equations or systems of equations of order $$\alpha$$. All that is needed is the RHS of the equation and the Jacobian. Because your equation is formulated with a term in $$D^v$$ and one in $$D^{2\cdot v}$$, given that $$D^v D^v = D^{2\cdot v}$$, you can rewrite your single equations as the following system of two equations: $$\begin{equation*} \left\{ \begin{array}{lcr} D^v y(x)=z(x),\\ \alpha\cdot D^{v}z(x) + \beta\cdot z(x) + \gamma\cdot y(x) = f(x). \end{array} \right. \end{equation*}$$

I hope this helps.

• Simply a comment that it would be really interesting, at least to me, to see what solutions for the 'fractional simple harmonic oscillator' $\alpha=1$, $\beta=0$, $\gamma>0$, for example, look like. Apr 6, 2023 at 0:43
• Would you mind doing this for a specific example of this fde and adding it in your answer? $$\qquad$$ If you do that, I'll also give you the bounty (unless someone else comes along that offers a non-numeric solution). Because just like that, the answer seems too little and unspecific to me. $$\qquad$$ So I would be grateful for an example. Apr 8, 2023 at 0:58
• @Aruralreader, I don't think you'd have a FC harmonic oscillator for v<0.5. Apr 11, 2023 at 1:40
• @gciriani, I am interested and naive, so thank you. Apr 11, 2023 at 1:47

You can factor the equation into two operators,

$$(D^{2v}+aD^v+b)u=\left(D^v+\frac{-a+\sqrt{a^2-4b}}{2}\right)\left(D^v+\frac{-a-\sqrt{a^2-4b}}{2}\right)u,$$ assuming $$a,b$$ are constants since $$D^v$$ is linear and since $$D^{u+v}=D^uD^v.$$ Now you only need to solve $$D^vu-\alpha u=w,$$ for arbitrary constant $$\alpha.$$ The only issue now is... what is $$v?$$ If it's $$1/2$$, notice that

$$(D^{\frac{1}{2}}-\alpha)u=w\iff (D^{\frac{1}{2}}+\alpha)(D^{\frac{1}{2}}-\alpha)u=(D^{\frac{1}{2}}+\alpha)w\\ \iff (D-\alpha^2)u=D^{\frac{1}{2}}w+\alpha w.$$ And you can proceed from there, since $$w$$ is known. Now, if $$v=p/q$$, then simply apply $$\prod_{k=1}^{q-1}\left(D^{\frac{p}{q}}-\alpha\omega^k_q\right)$$ to both sides ($$\omega_q$$ being the $$q$$th roots of unity) to obtain $$\left(D^p-\alpha^q\right)u=w^*,$$ $$w^*$$ meaning $$\prod\left(D^{\frac{p}{q}}-\alpha\omega^k_q\right)w$$, and you are left with solving for $$u$$. Fortunately, you can also decompose $$D^p-\alpha^q=\prod_{k=0}^{p-1}\left(D-\omega_p^k\alpha^{q/p}\right),$$ and we obtain

$$\boxed{\boxed{u=\prod_{k=0}^{p-1}\left(D-\omega_p^k\alpha^{q/p}\right)^{-1}\prod_{k=1}^{q-1}\left(D^{\frac{p}{q}}-\alpha\omega^k_q\right)w}}$$ from $$(D^{\frac{p}{q}}-\alpha)u=w.$$ Just remember that $$(D+\alpha)^{-1}w=e^{-\alpha t}\left(\int e^{\alpha t}wdt+C_0\right)$$ and $$D^{\frac{1}{n}}u=\frac{1}{\left(-\frac{n+1}{n}\right)!}\int^t(t-\xi)^{-\frac{n+1}{n}}u(\xi)d\xi,$$ and you've got the solution. Also, $$\frac{1}{\left(-\frac{n+1}{n}\right)!}=\frac{\sin\left(\pi \frac{n+1}{n}\right)(\frac{1}{n})!}{\pi}$$ from $$\Gamma(z)\Gamma(1-z)=\frac{\pi}{\sin \pi z}.$$

Assuming that we accept the boundary condition $$\operatorname{D}^{\aleph}\left[ y\left( x \right) \right] \operatorname{\mid}\limits_{x \to 0^{-}} = y_{\aleph,\, 0}$$, then we could work with the definitions of Laplace transforms for fractional derivatives of a function. Let's say: \begin{align*} \mathcal{L}_{x}\left\{ \operatorname{D}^{\aleph}\left[ y\left( x \right) \right] \right\}\left( s \right) = s^{\aleph} \cdot \mathcal{L}_{x}\left\{ y\left( x \right) \right\}\left( s \right) - \sum\limits_{k = 1}^{\operatorname{round}\left( \left| \aleph \right| \right)}\left[ s^{\aleph - k} \cdot y_{k - 1,\, 0} \right] \wedge \aleph \in \mathbb{C} \setminus \left\{ \mathbb{Z} \setminus \mathbb{N} \right\}\\ \end{align*}

So: \begin{align*} f\left( x \right) =\, &\alpha \cdot \operatorname{D}^{2 \cdot v} \left( y\left( x \right) \right) + \beta \cdot \operatorname{D}^{v} \left( y\left( x \right) \right) + \gamma \cdot y\left( x \right)\\ \mathcal{L}_{x}\left\{ f\left( x \right) \right\}\left( s \right) =\, &\alpha \cdot \left( s^{2 \cdot v} \cdot \mathcal{L}_{x}\left\{ y\left( x \right) \right\}\left( s \right) - \sum\limits_{k = 1}^{\operatorname{round}\left( \left| 2 \cdot v \right| \right)}\left[ s^{2 \cdot v - k} \cdot y_{k - 1,\, 0} \right] \right) +\\ &\beta \cdot \left( s^{v} \cdot \mathcal{L}_{x}\left\{ y\left( x \right) \right\}\left( s \right) - \sum\limits_{k = 1}^{\operatorname{round}\left( \left| v \right| \right)}\left[ s^{v - k} \cdot y_{k - 1,\, 0} \right] \right) +\\ &\gamma \cdot \mathcal{L}_{x}\left\{ y\left( x \right) \right\}\left( s \right)\\ \mathcal{L}_{x}\left\{ y\left( x \right) \right\}\left( s \right) =\, &\frac{\left( \alpha - \beta \right) \cdot \left( -\sum\limits_{k = 1}^{\operatorname{round}\left( \left| 2 \cdot v \right| \right)}\left[ s^{2 \cdot v - k} \cdot y_{k - 1,\, 0} \right] - \sum\limits_{k = 1}^{\operatorname{round}\left( \left| v \right| \right)}\left[ s^{v - k} \cdot y_{k - 1,\, 0} \right] \right) - \mathcal{L}_{x}\left\{ f\left( x \right) \right\}\left( s \right)}{\alpha \cdot s^{2 \cdot v} + \beta \cdot s^{v} + \gamma}\\ y\left( x \right) =\, &\mathcal{L}_{s}^{-1}\left\{ \frac{\left( \alpha + \beta \right) \cdot \left( -\sum\limits_{k = 1}^{\operatorname{round}\left( \left| 2 \cdot v \right| \right)}\left[ s^{2 \cdot v - k} \cdot y_{k - 1,\, 0} \right] - \sum\limits_{k = 1}^{\operatorname{round}\left( \left| v \right| \right)}\left[ s^{v - k} \cdot y_{k - 1,\, 0} \right] \right) - \mathcal{L}_{x}\left\{ f\left( x \right) \right\}\left( s \right)}{\alpha \cdot s^{2 \cdot v} + \beta \cdot s^{v} + \gamma} \right\}\left( x \right)\\ \end{align*}

So the solution would be: \fbox{ \begin{align*} y\left( x \right) &= -\mathcal{L}_{s}^{-1}\left\{ \frac{\left( \alpha + \beta \right) \cdot \left( \sum\limits_{k = 1}^{\operatorname{round}\left( \left| 2 \cdot v \right| \right)}\left[ s^{2 \cdot v - k} \cdot y_{k - 1,\, 0} \right] + \sum\limits_{k = 1}^{\operatorname{round}\left( \left| v \right| \right)}\left[ s^{v - k} \cdot y_{k - 1,\, 0} \right] \right) + \mathcal{L}_{x}\left\{ f\left( x \right) \right\}\left( s \right)}{\alpha \cdot s^{2 \cdot v} + \beta \cdot s^{v} + \gamma} \right\}\left( x \right) \end{align*} }

Of course, the solution has a problem, besides being not a closed form: the definition of the laplace transform is Libra, so the solution only applies to a certain fractional derivative operator but it's still a sloution.