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

Question

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.

Answer 1

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.

Answer 2

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}.$

Answer 3

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.