Метод Рунге-Кутта

[Maik2013]

Здравствуйте, пожалуйста помогите в программе Delphi решить систему обыкновенных дифференциальных уравнений методом Рунге-Кутта и вывести график. Для одного уравнения нашёл код, для моего уравнения как выглядит код программа.

$$ \begin{equation} \left\{ \begin{array}{ll} \dfrac{dx}{dt}=-2x+4y,\\ \dfrac{dy}{dt}=-x+3y \end{array} \right. \end{equation} $$

начальными условиями $x(0)=3,$ $y(0)=0.$

частное решение:

$$\begin{equation}\label{9} \left\{ \begin{array}{ll} x(t)=4e^{-t}-e^{2t},\\ y(t)=e^{-t}-e^{2t}.\\ \end{array} \right. \end{equation}$$

Код программа:

Используется синтаксис Delphi

unit Unit1;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, TeEngine, Series, ExtCtrls, TeeProcs, Chart;
 
type
  TForm1 = class(TForm)
    Chart1: TChart;
    Series1: TLineSeries;
    Memo1: TMemo;
    Edit1: TEdit;
    Edit2: TEdit;
    Edit3: TEdit;
    Edit4: TEdit;
    Button1: TButton;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Button2: TButton;
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
      private
    { Private declarations }
  public
    { Public declarations }
  end;
 
var
  Form1: TForm1;
    y0, x0, a, b, h: extended;
 
implementation
 
{$R *.dfm}
   function fxy(x,y: extended): extended;
begin
 fxy:=-(2*x+ y)/(1+cos(y));
end;
 
procedure TForm1.Button1Click(Sender: TObject);
var
 x,y : extended;
 k1,x1,y1,k2,y2,
   k3,y3,x2,k4,dy       :extended;
begin
y0 := StrToFloat(Edit1.Text);
x0 := 0;
 a  := StrToFloat(Edit2.Text);
 b  := StrToFloat(Edit3.Text);
 h  := StrToFloat(Edit4.Text);
 memo1.Lines.Add('Метод Рунге-Кутта');
 x:=x0;
 y:=y0;
 memo1.Lines.Add('  x='+floattostr(x)+';                y='+floattostr(y));
 Chart1.SeriesList[0].AddXY(x,y,'',clblue);
 while x<=b do
  begin
 
   k1:=h*fxy(x,y);
   x1:=x+h/2;
   y1:=y+k1/2;
 
   k2:=h*fxy(x1,y1);
   y2:=y+k2/2;
 
   k3:=h*fxy(x1,y2);
   y3:=y+k3;
 
        x2:=x+h;
   k4:=h*fxy(x2,y3);
   dy:=(k1+2*k2+2*k3+k4)/6;
   y:=y+dy;
   x:=x+h;
   memo1.Lines.Add('  x='+floattostr(x)+';              y='+floattostr(y));
   Chart1.SeriesList[0].AddXY(x,y,'',clblack);
  end;
 

Еще, пропуская несколько значение $x$, $y$ потом выходил результат на Memo какой команда можно использовать?

Помогите пожалуйста, очень нужны.

[Pavia]

Это не уравнение. Это система уравнений. Для её решения вам надо записать её в нормальной форме Коши. Т.е. с использованием обозначений $x_1$, $x_2$, $f_1(t,x_1,x_2,...)$, $f_2(t,x_1,x_2,...)$
На каждом шаге интегрирования применять расчет нового значения каждой функции как в методе Рунге-Кутты используя опорные переменные(внутри шага отдельно к каждому уравнению). В конце шага присваиваешь полученные результаты, опорным переменным для всей системы.

Цитата:

Еще, пропуская несколько значение , потом выходил результат на Memo какой команда можно использовать?

Используйте условный оператор "if". К примеру:

 if not((5<x) and (x<=10)) then memo1.Lines.Add('  x='+floattostr(x)+';              y='+floattostr(y)); 

Пропускаем все что между 5 и 10.

[Maik2013]

Pavia

Pavia в сообщении #1026929 писал(а):

if not((5<x) and (x<=10)) then memo1.Lines.Add(' x='+floattostr(x)+'; y='+floattostr(y));

Спасибо за это.


Вы имейте виду так

$$ \dfrac{dx}{dt}=-\dfrac{d^2y}{dt^2}+3\dfrac{dy}{dt} $$

[nazarov_m]

А вам обязательно нужно решать это в Дельфи? Просто из наследников паскаля актуальным на данный момент считается одна Ада.
Ваша задача на Аде вполне решалась и в общем векторном случае:

(Оффтоп)

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Numerics.Generic_Real_Arrays;

procedure RungeKuttaDemo is

   type Real is digits 15;
   package FIO is new Ada.Text_IO.Float_IO(Real); use FIO;
   
   package Float_arrays is new Ada.Numerics.Generic_Real_Arrays(Real); 
   use Float_arrays;
   package Math is new Ada.Numerics.Generic_Elementary_Functions(Real);
   
   type Derivative is access function( t:  Real; y:  Real_Vector)  return Real_Vector;
   
   ----------------------------------------------------------
   
   function Vector(A:   Real_Matrix;  i: Integer) return Real_Vector is
      V:    Real_Vector(A'Range(1));
   begin
      for j in A'Range(1) loop
         V(j) := A(j,i);
      end loop;
      return V;
   end Vector;
   
   function CoVector(A:  Real_Matrix; j: Integer) return Real_Vector is
      V:    Real_Vector(A'Range(2));
   begin
      for i in A'Range(2) loop
         V(i) := A(j,i);
      end loop;
      return V;
   end CoVector;
   
   Pragma Inline(Vector);
   Pragma Inline(CoVector);
   
   -----------------------------------------------------------
   
   procedure Runge (   yp_func : Derivative; 
               t :  in out Real_Vector;
               y_value: in out  Real_Matrix;
               dt : Real) is
      ------------------
      dy1, dy2, dy3, dy4 : Real_Vector(y_value'Range(1));
      y: Real_Vector(y_value'Range(1));
      -------------------
   begin
      for n in t'First .. t'Last-1 loop--------------
      y := Vector(y_value,n);
      ------------
      dy1 := dt * yp_func(t(n), y);
      dy2 := dt * yp_func(t(n) + dt / 2.0, y + dy1 / 2.0);
      dy3 := dt * yp_func(t(n) + dt / 2.0, y + dy2 / 2.0);
      dy4 := dt * yp_func(t(n) + dt, y + dy3);
      t(n+1) := t(n) + dt;
      ------------    
      for l in y_value'Range(1) loop   
         y_value(l,n+1) := y(l) + 
            (dy1(l) + 
               2.0 * (dy2(l) + dy3(l)) + dy4(l)) / 6.0;
      end loop;
      end loop;--------------------------------------
   end Runge;
   
   ----------------------------------------------------------------
 
   procedure Print (Name: Character; t, y : Real_Vector; modnum : Positive) is 
   begin
      for i in t'Range loop
         if i mod modnum = 0 then
         Put(Name&"(");   Put (t(i), Exp=>0, Fore=>0, Aft=>1);
         Put(") = "); Put (y(i), Exp=>0, Fore=>0, Aft=>8);
         New_Line;
         end if;
      end loop;
   end Print;
   
   ----------------------------------------------------------------
   
   function ODE_System_Example   (t:  Real; y:  Real_Vector )  return Real_Vector is
      y_new : Real_Vector(y'Range);
   begin
      y_new(1) := -2.0*y(1)    +    4.0*y(2);
      y_new(2) := -y(1)      +   3.0*y(2);   
      -------
      return y_new;
   end ;
   
   ---------------------------------------------------------------

   dt : constant Real := 0.10;
   N : constant Positive := 100;
   t_arr : aliased Real_Vector(0 .. N);
   y_matrix: aliased Real_Matrix(1 .. 2, 0 .. N );
   ---- y(1) -- will be your x 
   ---- y(2) -- will be your y
   
begin
   t_arr(0) := 0.0;
   
   y_matrix(1,0) := 3.0;  --- x0
   y_matrix(2,0) := 0.0;  --- y0

   Runge (   ODE_System_Example'Access,
         t_arr,
         y_matrix,
         dt);
         
   Print ('x',t_arr, CoVector(y_matrix,1), 10);
   Print ('y',t_arr, CoVector(y_matrix,2), 10);
   
end RungeKuttaDemo;

[Maik2013]

nazarov_m

Да потому, что чуть чуть разбираюсь.

Я другой задачу решу, проста это для контроля, программа работает или нет.

[nazarov_m]

Maik2013 в сообщении #1028790 писал(а):

nazarov_m

Да потому, что чуть чуть разбираюсь.

Я другой задачу решу, проста это для контроля, программа работает или нет.

Понятно.