1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Solving a fractional differential equation numerically

  1. Mar 7, 2013 #1
    1. The problem statement, all variables and given/known data

    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 ([itex]\varepsilon[/itex] below) defined for each of my time points, t. I calculate the fractional derivative of the strain ([itex]D^{\alpha}\varepsilon[/itex] below) using a numerical method according to the Grünwald-Letnikov definition.

    I want to solve for the stress ([itex]\sigma[/itex] below).

    2. Relevant equations

    The differential equation that describes the relationshp between strain ([itex]\varepsilon[/itex]) and stress ([itex]\sigma[/itex]) is:


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

    3. The attempt at a solution

    My naive approach was to rearrange the equation to yield:


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

    [tex]D^{\alpha}\sigma\left(t\right)=E\left(D^{\alpha} \varepsilon \left(t\right)-\frac{\sigma\left(t-1\right)}{\eta}\right)[/tex]

    And integrate with:


    I suspect there may be two problems with this approach, but I am not certain:
    • First, the calculation of the fractional derivative uses [itex]\sigma\left(t-1\right)[/itex] instead of [itex]\sigma\left(t\right)[/itex]. 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 [itex]D^{\alpha}\sigma[/itex] to yield the resulting stress [itex]\sigma[/itex]? 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.
  2. jcsd
  3. Mar 7, 2013 #2
    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):

    Code (Text):

    % 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);
    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 [Broken]
    Last edited by a moderator: May 6, 2017
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted