Уравнения в частных производных

        function norm2psi: real;

        function residual1: real;

        function residual2: real;

        function a(k: integer): real;

        function b(k: integer): real;

        procedure Button1Click(Sender: TObject);

      private

        { Private declarations }

      public

        { Public declarations }

      end;

var

      Form1: TForm1;

      N: integer;

      m: integer;

      x: array of real;

      alfa: array of real;

      beta: array of real;

implementation

{$R *.dfm}

procedure TForm1.FormCreate(Sender: TObject);

    begin

    ComboBox1.Items.Add('~0.2'); 

    ComboBox1.Items.Add('~0.01');

    ComboBox1.Items.Add('~0.001');

    ComboBox1.ItemIndex:=0;

    RadioGroup1.Items.Add('phi(x)=exp(x)-1-exp(Pi/2)x; psi(x)=x*x');

    RadioGroup1.Items.Add('phi(x)=ln(x+1)-2/(Pi+2)x; psi(x)=exp(x)-1');

    RadioGroup1.ItemIndex:=0;

    RadioGroup2.Items.Add('Метод средних');

    RadioGroup2.Items.Add('Модифицированный  метод средних');

    RadioGroup2.ItemIndex:=0;

end;

    // граничные условия

    function TForm1.phi(t: real): real;

    begin

      case RadioGroup1.ItemIndex of

          0: result:=exp(t)-1-exp(Pi/2)*t;

          1: result:=ln(t+1)-2/(Pi+2)*t;

end;

end;

        function TForm1.psi(t: real): real;

    begin

      case RadioGroup1.ItemIndex of

          0: result:=t*t;

          1: result:=exp(t)-1;

 end;

end;

    // приближенное решение и ее производная по времени

    function TForm1.v(t: real): real;

    var i: integer;

    begin 

    result:=alfa[0]*sin(t);

    if m>1 then

    for i:=1 to m-1 do result:=result+alfa[i]*sin((2*i+1)*t);

    end;

function TForm1.vt(t: real): real;

    var i: integer;

    begin

    result:=beta[0]*sin(t);

    if m>1 then

    for i:=1 to m-1 do result:=result+(2*i+1)*beta[i]*sin((2*i+1)*t);

    end;

    //норма  в L2(a;b) функции phi(x)

function TForm1.norm1phi: real;

    var i: integer;

    begin

    result:=phi(x[0])*phi(x[0]);

    for i:= 1 to N-2 do result:=result+phi(x[i])*phi(x[i]);

    result:=sqrt(result);

    end;

    //норма в С функции phi(x)

function TForm1.norm2phi: real;

    var i: integer;

    begin

    result:=abs(phi(x[0]));

    for i:=1 to N-2 do

      if abs(phi(x[i]))>result then result:=abs(phi(x[i]));

    end;

    //норма  в L2(a;b) функции psi(x)

function TForm1.norm1psi: real; 

    var i: integer;

    begin

    result:=psi(x[0])*psi(x[0]);

    for i:= 1 to N-2 do result:=result+psi(x[i])*psi(x[i]);

    result:=sqrt(result);

    end;

    //норма в С функции psi(x)

function TForm1.norm2psi: real;

    var i: integer;

    begin

    result:=abs(psi(x[0]));

    for i:=1 to N-2 do

      if abs(psi(x[i]))>result then result:=abs(psi(x[i]));

    end;

    //-------------------------------------------------------------

function TForm1.r1phi: real;

    var i: integer;

    begin

    result:=Power(phi(x[0])-v(x[0]), 2);

    for i:=1 to N-2 do result:=result+Power(phi(x[i])-v(x[i]), 2);

    result:=sqrt(result);

    end;

function TForm1.r1psi: real;

    var i: integer;

    begin

    result:=Power(psi(x[0])-vt(x[0]), 2);

    for i:=1 to N-2 do result:=result+Power(psi(x[i])-vt(x[i]), 2);

    result:=sqrt(result); 
 

    end;

    //-----------------------------------------------------------------

function TForm1.r2phi: real;

    var i: integer;

    begin

    result:=abs(phi(x[0])-v(x[0]));

    for i:=1 to N-2 do

      if abs(phi(x[i])-v(x[i]))>result

       then result:=abs(phi(x[i])-v(x[i]));

end; 

function TForm1.r2psi: real;

    var i: integer;

    begin

    result:=abs(psi(x[0])-vt(x[0]));

    for i:=1 to N-2 do

      if abs(psi(x[i])-vt(x[i]))>result

       then result:=abs(psi(x[i])-vt(x[i]));

end;

    //невязка 1 норма в L

    function TForm1.residual1: real;

    begin

    result:=(r1phi/norm1phi+r1psi/norm1psi)*100;

    end;

    //невязка  2 (норма в C[a;b])

    function TForm1.residual2: real;

    begin

    result:=(r2phi/norm2phi+r2psi/norm2psi)*100;

end; 

    //вычисление коэф-ов

function TForm1.a(k: integer): real;

    var i: integer;

    begin

      case RadioGroup2.ItemIndex of

        0: begin

            result:=0;

            for i:=0 to N-2 do

             result:=result+phi(x[i])*sin((2*k+1)*x[i])*Pi/2/(N-1);

             result:=result*4/Pi;

end;

        1: begin

            result:=0;

            for i:=0 to N-2 do

    result:=result+phi(x[i])*(cos((2*k+1)*(i+1)*Pi/2/(N-1))-cos((2*k+1)*i*Pi/2/(N-1)));

    result:=result*(-4)/(2*k+1)/Pi;

            end;

      end;

    end; 

function TForm1.b(k: integer): real;

    var i: integer;

    begin

      case RadioGroup2.ItemIndex of

        0: begin

            result:=0;

            for i:=0 to N-2 do

             result:=result+psi(x[i])*sin((2*k+1)*x[i])*Pi/2/(N-1);

             result:=result*4/(2*k+1)/Pi;

           

     end;

        1: begin

            result:=0;

            for i:=0 to N-2 do

    result:=result+psi(x[i])*(cos((2*k+1)*(i+1)*Pi/2/(N-1))-cos((2*k+1)*i*Pi/2/(N-1)));

    result:=result*(-4)/(2*k+1)/(2*k+1)/Pi;

            end;

      end;

    end;

procedure TForm1.Button1Click(Sender: TObject);

    var i: integer;

    begin

    // количество узлов разбиения промежутка  в зависимости от шага

      case ComboBox1.ItemIndex of

        0: N:=8;

        1: N:=157;

        2: N:=1571;

      end;

    //заполняем  массив x[i]

    SetLength(x, N);

    for i:=0 to N-1 do x[i]:=Pi/2/(N-1)*i;

    for i:=0 to N-2 do x[i]:=(x[i+1]+x[i])/2;

    //заполняем  массив alfa[i] и beta[i]

    SetLength(alfa, 101);

    for i:=0 to 100 do alfa[i]:=a(i); 

    SetLength(beta, 101);

    for i:=0 to 100 do beta[i]:=b(i);

    // на график 

    Chart1.Series[0].Clear;

    Chart1.Series[1].Clear; 

    for i:=1 to 100 do begin  //кол-во слагаемых

                        m:=i;

                        Chart1.Series[0].Add(residual1, IntToStr(m));

                        Chart1.Series[1].Add(residual2, IntToStr(m));

                        end;

    {m:=70;

    for i:=0 to 157 do begin

                        Chart1.Series[0].Add(psi(i/100));

                        Chart1.Series[1].Add(vt(i/100));

                        end;} 

    for i:=0 to 10 do Memo1.Lines.Add('');

    for i:=0 to 10 do Memo1.Lines[i]:='a'+IntToStr(i)+'='+FloatToStr(RoundTo(a(i), -5)); 

    for i:=0 to 10 do Memo2.Lines.Add('');

    for i:=0 to 10 do Memo2.Lines[i]:='b'+IntToStr(i)+'='+FloatToStr(RoundTo(b(i), -5)); 

    Edit1.Text:=FloatToStr(norm1phi);

    Edit2.Text:=FloatToStr(norm1psi);

    Edit3.Text:=FloatToStr(norm2phi);

    Edit4.Text:=FloatToStr(norm2psi);

    Edit5.Text:=FloatToStr(residual1);

    Edit6.Text:=FloatToStr(residual2);

    end;

    end.