Здравствуйте. Нужен взгляд человека понимающего 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):
Ik+1=Ik+Δt*(U-Ik-R-ωk*Ce)/L
ωk+1k+Δt*(Ik*Cm-Mст)/J

и собственно сам метод Рунге-Кутты
k11 = dt*((U-Ik*R-ωk*Ce)/L)
k21 = dt*((U-Ik*R-ωk*Ce)/L) + k11/2
k31 = dt*((U-Ik*R-ωk*Ce)/L) + k21/2
k41 = dt*((U-Ik*R-ωk*Ce)/L) + k31/2

k12 = dt*((Ik*Cm-Mc)/J)
k22 = dt*((Ik*Cm-Mc)/J) + k12/2
k32 = dt*((Ik*Cm-Mc)/J) + k22/2
k42 = dt*((Ik*Cm-Mc)/J) + k32/2
тогда
Ik=(Ik+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.