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