Приближенное вычисление левосторонней дробной производной Римана-Лиувилля

      Введем  обозначения

 

      Наиболее  простой алгоритм образуется из определения  Грюнвальда-Летникова 

А именно, опустив , получим 

С помощью  рекурсивной формулы  

алгоритм (5) при может быть осуществлен по схеме «умножение-сложение-…»: 

который позволяет избежать прямого использования гамма-функции. 
 
 
 
 

Вычисление  теплового потока. 

     С точки зрения практики и технологии, одним из важнейших контролируемых параметров печи с форсированной  тягой является интенсивность теплового  потока в ее стенках (см. [1]). Стандартный метод решения этой задачи состоит в вычислении теплового потока посредством измерения двух температур в стенке, последовательное воспроизведение температурного поля и вычисление теплового потока в некоторой точке через разность температур в двух соседних точках. Этот метод требует два измерения температуры в двух точках, находящихся на различной глубине стенки печи. Мы рассмотрим другой подход, основанный на использовании дробных производных Грюнвальда-Летникова.

     Рассмотрим  следующее одномерное уравнение  теплопроводности для полубесконечного тела

     , t > 0                     (6) 
 
 

     где – время, пространственная координата направления теплопроводности, – теплоемкость, - массовая плотность, - температура, – коэффициент теплопроводности. Запишем вспомогательную функцию 

     которая является решением задачи

       (7) 
 
 

     Преобразование  Лапласа уравнения (6) дает 

     Решение уравнения (7), ограниченное при имеет вид 

     откуда  мы найдем 

     Получим 

     Применяя  обратное преобразование Лапласа получим 

 и  используя свойство () и линейность дробного интегродифференцирования мы придем к выражению 

Теперь  мы можем вернуться к функции с помощью соотношения Принимая во внимание условие мы получим основное аналитическое отношение для вычисления теплового потока в точке

   где , а

 искомый  тепловой поток. 

Рассмотрим  функцию  

Точное  значение производной  
 

Текст программы

unit Unit1; 

interface 

uses

  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

  Dialogs, Gamma, StdCtrls, Grids, TeEngine, Series, ExtCtrls, TeeProcs, Chart; 
 

type

  TForm1 = class(TForm)

    StringGrid1: TStringGrid;

    Button1: TButton;

    Label1: TLabel;

    Label2: TLabel;

    Label3: TLabel;

    Label4: TLabel;

    Label5: TLabel;

    Edit1: TEdit;

    Edit2: TEdit;

    Edit3: TEdit;

    Edit4: TEdit;

    Edit5: TEdit;

    Chart1: TChart;

    Series1: TPointSeries;

    Series2: TPointSeries;

    procedure Button1Click(Sender: TObject);

   

  private

    { Private declarations }

  public

    { Public declarations }

  end; 

var

  Form1: TForm1; 

implementation 

{$R *.dfm}

function f(x:real):real;

begin

  f:=x;

end; 
 
 

procedure TForm1.Button1Click(Sender: TObject);

  const n_n=100;    {определяет с какой точностью нах-ся реш.}

  var i1,i2,k,j,i:integer;

      x,h,sum:real;

      toch,pribl,Ea,Eo:real;

      nom:integer;

      w:array[0..n_n] of real;

      al,b,d,p:real;

      n:integer; 

begin

  n:=strtoint(edit1.text);

  al:=strtofloat(edit2.text);

  b:=strtofloat(edit3.text);

  d:=strtofloat(edit4.text);

  p:=strtofloat(edit5.text);

  h:=(d-b)/n;

  stringgrid1.RowCount:=1;

  stringgrid1.ColCount:=7;

  stringgrid1.Cells[0,0]:='n';

  stringgrid1.Cells[1,0]:='x';

  stringgrid1.Cells[2,0]:='f(x)';

  stringgrid1.Cells[3,0]:='точное';

  stringgrid1.Cells[4,0]:='приближенное';

  stringgrid1.Cells[5,0]:='абс.погр.';

  stringgrid1.Cells[6,0]:='относит.погр.';

  x:=b;

  k:=0;

  j:=1;

  nom:=0;

  While abs(x-d-h)>0.000001 do

    begin

      stringgrid1.RowCount:=stringgrid1.RowCount+1;

      k:=k+1;

      sum:=0;

      toch:=(GammaF2(p+1)*exp((p-al)*ln(x)))/GammaF2(p-al+1);

      w[0]:=1;

      For i:=1 to n_n do w[i]:=(1-(al+1)/i)*w[i-1];

      For i:=0 to n_n do sum:=sum+w[i]*f(x-i*(x/n_n));

      pribl:=sum*exp((-al)*ln(n_n/x));

      series1.AddXY(x,toch);

      series2.AddXY(x,pribl);

      Ea:=abs(pribl-toch);

      Eo:=Ea/toch;

      stringgrid1.Cells[0,j]:=inttostr(nom);

      stringgrid1.Cells[1,j]:=floattostr(x);

      stringgrid1.Cells[2,j]:=floattostr(f(x));

      stringgrid1.Cells[3,j]:=floattostr(toch);

      stringgrid1.Cells[4,j]:=floattostr(pribl);

      stringgrid1.Cells[5,j]:=floattostr(Ea);

      stringgrid1.Cells[6,j]:=floattostr(Eo);

      nom:=nom+1;

      x:=x+h;

      k:=k+1;

      j:=j+1;

    end; 
 
 
 

end; 

end. 

unit Gamma;

interface

uses Math, Sysutils; 

   function GammaF2(x : Double):Double;

   function GammaObr(x : Double):Double;

implementation 

   function GammaF2(x : Double):Double;

   var G:Double;

           GamF: Double;

   const        A1=-0.422784335092;

                A2=-0.233093736365;

                A3=0.191091101162;

                G1=-0.00000018122;

                G2=0.000001328554;

                G3=-0.000002625721;

                G4=-0.000017527917;

                G5=0.000145624324;

                G6=-0.000360851496;

                G7=-0.000804341335;

                G8=0.008023278113;

                G9=-0.017645242118;

                G10=-0.024552490887;

           begin 

begin

          GamF:=1;

    WHILE  x>1 DO

      BEGIN

            x:=x-1;

        GamF:=GamF*x;

      END;

    GamF:=1/GamF;

      WHILE X<-1 DO

      BEGIN

            GamF:=GamF*X;

            X:=X+1;

      END;

    IF X<>1 THEN

      BEGIN

     G:=(((((((((G1*X+G2)*X

             +G3)*X+G4)*X+

            G5)*X+G6)*X+G7)*X+G8)*X

          +G9)*X+G10)*X; 

 GamF:=GamF*((((G+A3)*X+A2)*X+A1)*X+1)*X*(1+X); 

  END;

          GamF:=1/GamF;

      GammaF2:=GamF;

  end;

end; 

function GammaObr(x : Double):Double;

   var G:Double;

           GamF: Double;

   const        A1=-0.422784335092;

                A2=-0.233093736365;

                A3=0.191091101162;

                G1=-0.00000018122;

                G2=0.000001328554;

                G3=-0.000002625721;

                G4=-0.000017527917;

                G5=0.000145624324;

                G6=-0.000360851496;

                G7=-0.000804341335;

                G8=0.008023278113;

                G9=-0.017645242118;

                G10=-0.024552490887;

           begin 

begin

          GamF:=1;

    WHILE  x>1 DO

      BEGIN

            x:=x-1;

        GamF:=GamF*x;

      END;

    GamF:=1/GamF;

      WHILE X<-1 DO

      BEGIN

            GamF:=GamF*X;

            X:=X+1;

      END;

    IF X<>1 THEN

      BEGIN

     G:=(((((((((G1*X+G2)*X

             +G3)*X+G4)*X+

            G5)*X+G6)*X+G7)*X+G8)*X

          +G9)*X+G10)*X; 

 GamF:=GamF*((((G+A3)*X+A2)*X+A1)*X+1)*X*(1+X); 

  END;

             GammaObr:=GamF;

  end;

end;

end. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Литература.

1. Podlubny I. Fractional Differential Equatious.-NewYork-London, Academic Press, 1999, 365 p.
2. Oldham K.B., Spanier J. The fractional calculus-NewYork: Academic Press,1974

Academic Press, Inc. 1250 Sixtu Avenue, Sun Diegu, California, 1974, 240 p.

3. С.Г. Самко, А.А. Килбас, О.И. Маричев. Интегралы и производные дробного порядка и некоторые их приложения. 1987
4. Ландау Л.Д., Лифшиц Е.М. Механика сплошных сред.-М.: ГИТТЛ, 1954 100с.