Метод Рунге-Кутты 4го порядка, Метод Рунге-Кутты 4го порядка точности для решения уравнений ДПТ по вт Опции

[justlikenirvanasback]

Здравствуйте. Нужен взгляд человека понимающего Pascal. Расчёт уравнений двигателя постоянного тока по второму закону Киргофа методом Рунге-Кутты 4го порядка, кратко о нём:
2й Закон Киргофа для ДПТ:
U=I*R+L*∂I/∂t+ω*Ce;
I*Cm=J*∂ω/∂t+Mст
далее уравнения приводятся к форме задачи Коши:
∂I/∂t=(U-I-R-ω*Ce)/L
∂ω/∂t=(I*Cm-Mст)/J
далее уравнения дифференцируют для удобства записи (на шаге k+1):
I_k+1=I_k+Δt*(U-I_k-R-ω_k*Ce)/L
ω_k+1=ω_k+Δt*(I_k*Cm-Mст)/J
и собственно сам метод Рунге-Кутты
k11 = dt*((U-I_k*R-ω_k*Ce)/L)
k21 = dt*((U-I_k*R-ω_k*Ce)/L) + k11/2
k31 = dt*((U-I_k*R-ω_k*Ce)/L) + k21/2
k41 = dt*((U-I_k*R-ω_k*Ce)/L) + k31/2

k12 = dt*((I_k*Cm-Mc)/J)
k22 = dt*((I_k*Cm-Mc)/J) + k12/2
k32 = dt*((I_k*Cm-Mc)/J) + k22/2
k42 = dt*((I_k*Cm-Mc)/J) + k32/2
тогда
I_k=(I_k+1/6*(k11+k21*2+k31*2+k41))
ω_k=(ω_k+1/6*(k12+k22*2+k32*2+k42))


Теперь пытаюсь осуществить этот метод на Pascal. В 20 строке вещественное переполнение, могут быть и другие ошибки, поэтому прошу вашей помощи.

Program RK_4;

uses crt;

var i,n:integer;

    Ua,Ra,La,Cef,Cmf,Js,Mc:real;

    fi_,fw_:array[0..1000] of real;

    k1,k2,k3,k4,k11,k21,k31,k41,dt,t,tmax:real;

    ff:text;


function difurI(Ua,fi_,Ra,Cef,fw_,La:real):real;

begin

difurI:=(Ua-fi_*Ra-Cef*fw_)/La;

end;

function difurW(Cmf,fi_,Mc,Js:real):real;

begin

difurW:=(Cmf*fi_-Mc)/Js;

end;

begin

Ua:= 220.0;

Ra:= 0.03;

Cef:= 0.2;

Cmf:= 1.2;

La:= 0.01;

Js:= 1.0;

Mc:= 1000.0;


fi_[0]:=0;

fw_[0]:=0;

t:=0;

dt:=0.001;

n:=0;

write('Calculating Runge-Kutta4...');

while (n<100) do begin

k1:=dt*difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La);

k2:=dt*(difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La)+k1/2);

k3:=dt*(difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La)+k2/2);

k4:=dt*(difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La)+k3);

k11:=dt*difurW(t+dt,fi_[n],Mc,Js);

k21:=dt*(difurW(t+dt,fi_[n],Mc,Js)+k11/2);

k31:=dt*(difurW(t+dt,fi_[n],Mc,Js)+k21/2);

k41:=dt*(difurW(t+dt,fi_[n],Mc,Js)+k31);

t:=t+dt; inc(n);

fi_[n]:=fi_[n-1]+(k1+k2*2+k3*2+k4)/6;

fw_[n]:=fw_[n-1]+(k11+k21*2+k31*2+k41)/6;

end;writeln('Done!');


begin
clrscr;
assign(ff,'n.txt');
rewrite(ff);

      tmax:=100.0;
      dt:=0.001;

      for i:=1 to n do
      begin
        fi_[0]:=0.0;
      END;
      fi_[1]:=999;


      T:=0.0;
      writeln(ff,'  (t)        (i)       (w);');


while (t <= tmax) do
      begin
        difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La);
      end;
      begin
        difurW(t+dt,fi_[n],Mc,Js);
      end;
        T:=T+dt;
        write(ff,T:12:5);
        for i:=1 to n do
        begin
          write(ff,fi_[n]:12:5);
        end;
        begin
          write(ff,fw_[n]:12:5);
        end;

        writeln(ff,'');
      end;

close(ff);

END.