Solving a fractional differential equation numerically

[ajptech]

Homework Statement

This is a research problem that I imagine is very similar to a homework problem. I am a PhD student in biology, and I lack the mathematical background needed to make sense of the topic.

I would like to find the solution to the fractional differential equation that describes the behavior of a viscoelastic solid. The behavior can be described by a spring in series with a fractional-order dashpot (a "spring-pot").

I have the input strain ($\varepsilon$ below) defined for each of my time points, t. I calculate the fractional derivative of the strain ($D^{\alpha}\varepsilon$ below) using a numerical method according to the Grünwald-Letnikov definition.

I want to solve for the stress ($\sigma$ below).

Homework Equations

The differential equation that describes the relationshp between strain ($\varepsilon$) and stress ($\sigma$) is:

$$D^{\alpha}\varepsilon=\frac{\sigma}{\eta}+\frac{1}{E}D^{\alpha}\sigma$$

Where

• $\sigma$ is the viscoelastic stress, which is what I want to solve for.
• $\alpha$ is the order of the differentiation (I have this value).
• $\varepsilon$ is the viscoelastic strain (I have these value).
• $\eta$ is the viscosity of the spring-pot (I have this value).
• $E$ is the spring constant (I have this value).

The Attempt at a Solution

My naive approach was to rearrange the equation to yield:

$$D^{\alpha}\sigma=E\left(D^{\alpha}\varepsilon-\frac{\sigma}{\eta}\right)$$

Then, for each time point, I calculate the change in the stress by:

$$D^{\alpha}\sigma\left(t\right)=E\left(D^{\alpha} \varepsilon \left(t\right)-\frac{\sigma\left(t-1\right)}{\eta}\right)$$

And integrate with:

$$\sigma\left(t\right)=\sigma\left(t-1\right)+D^{\alpha}\sigma\left(t\right)$$

I suspect there may be two problems with this approach, but I am not certain:

• First, the calculation of the fractional derivative uses $\sigma\left(t-1\right)$ instead of $\sigma\left(t\right)$. Is this a valid approach for a numerical approximation, or is this a total no-no? If this is invalid, how can I approach the problem differently?
• Second, does it make sense to integrate the fractional derivative $D^{\alpha}\sigma$ to yield the resulting stress $\sigma$? I am not sure if this is a correct way to interpret a fractional derivative.

I would really appreciate some assistance! If it would make more sense to solve the problem analytically, that would suit me just fine... I just don't have the slightest clue how to do this, unfortunately.

 

[ajptech]

More details...

To provide just a bit more detail about what I am attempting to do, here is the MATLAB code I am using (excluding the definitions of the variables for the time vector, strain vector, and constants):

% Preallocate result vector
stress = zeros(size(t));

% Calculate fractional derivative of strain
DStrain = fderiv(alpha, strain, dt);

% Calculate stress step-by-step
for i = 2:numel(t);
    DStress(i) = E*(DStrain(i)-stress(i-1)/eta);
    stress(i) = stress(i-1) + DStress(i);
end

In the above code, fderiv() is a function found on MATLAB File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/13858-fractional-differentiator/content/fderiv.m