backward finite difference with fractional differential equat. (G-L definition)


by sedo_5
Tags: finite difference, fractional calculus, grunwald-letnikov, mathematica
sedo_5
sedo_5 is offline
#1
Sep21-11, 07:11 PM
P: 1
Hi everyone!
I'm having a problem with backward finite difference method. I'm dealing with this problem using Mathematica software (and I'm going to paste some code below, maybe you don't exactly know the syntax, but it's pretty easy I think)..

I need to use the backward finite difference method with the Grunwald-Letnikov fractional derivative definition (probably it's an uncommon topic..).
Long story short, G-L definition for the fractional derivative is this(two links here):


http://en.wikipedia.org/wiki/Gr%C3%B...kov_derivative
and
http://www.xuru.org/fc/GrunwaldLetnikov.asp

and can be usually used for backw.FD in order to approximate some fractional differential equation.

This is my problem.
I have a first order mass-spring differetial equation and and I have to write a code in Mathematica to approximate it (with backwFD method) and then compare it with the fractional differential equation case. (ps: the same problem appears in the second order equation (mass-spring-damper) approximation, but let's just limit this on the first order).

The problem is that, for the ordinary case I have no problem in using FD method, instead I have a big problem in using FD method for the fractional case: I'll write some code to be more precise.

---
first order - ordinary
---
massa = 0;cvisc = 5000;kappa = 1000;forzante[t_] := 2000*UnitStep[t];
LL = 30;
nn = 31;
deltah = LL/(nn - 1)
ii = Range[0, nn];
tt = deltah*(ii - 1);
Clear[u]; Clear[eq]; Clear[eqi]; Clear[ui];
ui = Array[u, nn + 1, 0];
eqi = Array[eq, nn + 1, 0];

For[i = 1, i <= nn, i++,
 eq[i] = cvisc*((u[i] - u[i - 1])/(deltah)) + kappa*(u[i]) == 
   forzante[tt[[i + 1]]]]
eq[0] = u[1] == 0;

solApprox := NSolve[eqi, ui]; uisol = 
 Table[ui[[i]] /. solApprox[[1]], {i, nn + 1}];
dataApprox = Table[{tt[[i]], uisol[[i]]}, {i, 2, nn + 1}];
plot1 = ListPlot[dataApprox, PlotStyle -> {Purple}, 
  PlotRange -> Automatic]

---
first order - fractional
---
everything as before, but with this FD cicle

For[i = 1, i <= nn, i++,
 eq[i] = (1/deltah)*cvisc*
     Sum[(-1)^r*Binomial[.5, r]*u[i - r], {r, 0, i}] + kappa*(u[i]) ==
    forzante[tt[[i + 1]]]]
which is made by substituting the first order derivative with a fractional derivative of order 1/2 using GL definition.

----

The first code (ordinary case) works flawlessly, but the second doesn't: the curve is completely different depending on the "deltah" imposed (it gives no problems when deltah=1 , but smaller values are generally used in approximations...and that's where this code gives problems).

Now, i've tried quite everything with this code but to no avail...the problem is for sure given by deltah but i cannot find a way.

It would be awesome if you could give me some hints on this. I think you could give me some hints also if you don't know Mathematica.

PS: Since it's my first post here (by the way...this forum is awesome and I say Hi to everybody :) ), please excuse me if I'm posting it in the wrong section or if I'm going against some netiquette.


Thanks a lot in advance.

R.
Phys.Org News Partner Science news on Phys.org
Going nuts? Turkey looks to pistachios to heat new eco-city
Space-tested fluid flow concept advances infectious disease diagnoses
SpaceX launches supplies to space station (Update)

Register to reply

Related Discussions
Finite Difference Discretization of a Fourth Order Partial Differential Term Differential Equations 4
Difference between Forward and Backward Fourier Transforms? Calculus 1
Numerical differentiation using forward, backward and central finite difference Differential Equations 2
System of ODE Boundary Value Problem with 2nd Order Backward Difference Differential Equations 1
Definition of Fractional Occupation Classical Physics 0