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;
Метод Рунге-Кутты 4го порядка Pascal
Здравствуйте. Нужен взгляд человека понимающего Pascal. Расчёт уравнений двигателя постоянного тока по второму закону Киргофа методом Рунге-Кутты 4го порядка, кратко о нём:
1.Закон Киргофа для ДПТ:
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+1=ωk+Δ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 строке вещественное переполнение, могут быть и другие ошибки, поэтому и прошу вашей помощи.
1.Закон Киргофа для ДПТ:
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+1=ωk+Δ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.
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.
Вместо типа Real использовать Double.
Цитата: ellor!
Вместо типа Real использовать Double.
Pascal ABC не знает double
Это зачем такое
Код:
for i:=1 to n do
begin
fi_[0]:=0.0;
END;
...
for i:=1 to n do
begin
write(ff,fi_[n]:12:5);
end;
begin
fi_[0]:=0.0;
END;
...
for i:=1 to n do
begin
write(ff,fi_[n]:12:5);
end;
Цитата: ellor!
Это зачем такое
Код:
for i:=1 to n do
begin
fi_[0]:=0.0;
END;
...
for i:=1 to n do
begin
write(ff,fi_[n]:12:5);
end;
begin
fi_[0]:=0.0;
END;
...
for i:=1 to n do
begin
write(ff,fi_[n]:12:5);
end;
что бы результаты вычислений записывались в файл
Цитата: Стасян Озёрный
Цитата: ellor!
Вместо типа Real использовать Double.
Pascal ABC не знает double
Попробуйте тип extended.
1) В основной программе La не инициализирована. Видимо по умолчанию инициализируется нулем или чем-то очень близким. При первом же вызове diffurI происходит деление на ноль.
2) Это что за ересь? Что это код, по вашему, должен делать?
2) Это что за ересь? Что это код, по вашему, должен делать?
Код:
странно,может я не там использую их
Код:
function difurI(Ua,fi_,Ra,Cef,fw_,La:real):double;
Цитата: P*t*
1) В основной программе La не инициализирована. Видимо по умолчанию инициализируется нулем или чем-то очень близким. При первом же вызове diffurI происходит деление на ноль.
2) Это что за ересь? Что это код, по вашему, должен делать?
2) Это что за ересь? Что это код, по вашему, должен делать?
1) как её инициализировать?
2)Код должен дифференцировать уравнения, пока время не достигнет максимального
Цитата: Стасян Озёрный
Цитата: P*t*
1) В основной программе La не инициализирована. Видимо по умолчанию инициализируется нулем или чем-то очень близким. При первом же вызове diffurI происходит деление на ноль.
2) Это что за ересь? Что это код, по вашему, должен делать?
2) Это что за ересь? Что это код, по вашему, должен делать?
1) как её инициализировать?
2)Код должен дифференцировать уравнения, пока время не достигнет максимального
1) Написать La := <число>
Это вы уж сами думайте, какое там число должно быть в соответствии с методом
2) Обратите внимание, difurI будет выполняться в цикле, но difurW находится уже после цикла и выполнится только один раз. Поизучайте синтаксис паскаля.
Код:
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;
begin
difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La);
end;
begin
difurW(t+dt,fi_[n],Mc,Js);
end;
1) Написать La := <число>
так там написано La:=0.01
2)В цикле это так?
[/quote]
или так?
так там написано La:=0.01
2)В цикле это так?
Код:
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;
begin
difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La);
end;
begin
difurW(t+dt,fi_[n],Mc,Js);
end;
или так?
Код:
while (t <= tmax) do
begin
difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La);
difurW(t+dt,fi_[n],Mc,Js);
end;
begin
difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La);
difurW(t+dt,fi_[n],Mc,Js);
end;
Цитата: Стасян Озёрный
1) Написать La := <число>
так там написано La:=0.01
так там написано La:=0.01
Ага, вот оно как... А зачем еще один begin перед clrscr ?
Цитата: Стасян Озёрный
2)В цикле это так?
Второй вариант правильный. А отступы в паскале вообще никакого значения не имеют, можно хоть всю программу в одну строчку написать.
можно и без него.
Всё равно вещественное переполнение в строке difurI:=(Ua-fi_*Ra-Cef*fw_)/La;
Код:
while (t <= tmax) do
begin
difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La);
difurW(t+dt,fi_[n],Mc,Js);
end;
begin
difurI(t+dt,fi_[n],fw_[n],Ua,Ra,La);
difurW(t+dt,fi_[n],Mc,Js);
end;
Это не может быть правильным, так как в этих двух строках не изменяются вообще никакие переменные и цикл будет бесконечным. Возвращаемые значения difurI и difurW также остаются неиспользованными.
end перед close(ff) вовсе не означает конец цикла, как вы могли от него ожидать, а является парным к begin перед clrscr
Ну, батенька, это полный бред сивой кобылы. Что вы написали в исходном коде.
Если отформатировать нормально код, то получится:
Сейчас нет времени разбираться с самим алгоритмом, но код ужасный и бессмысленный. Поэтому надо не искать
Если отформатировать нормально код, то получится:
Код:
Program RK_4;
uses crt;
var
i, n: integer;
Ua, Ra, La, Cef, Cmf, Js, Mc: real; //Cef, Cmf - не используются
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
// почему до 100, а не до 1000, если массивы объявлены так: fi_, fw_: array [0 .. 1000]
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 // зачем тут begin и соответственно внизу предпоследний end
clrscr; // сперва выводим на экран сообщения, и не давая прочитать пользователю, сразу стираем ???? - бред.
assign(ff, 'n.txt');
rewrite(ff);
tmax := 100.0;
dt := 0.001;
for i := 1 to n do
begin
fi_[0] := 0.0; // 100 раз присвоить fi_[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); //в цикл, если надо, не попало, + см выше по ф-ции difurI
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.
uses crt;
var
i, n: integer;
Ua, Ra, La, Cef, Cmf, Js, Mc: real; //Cef, Cmf - не используются
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
// почему до 100, а не до 1000, если массивы объявлены так: fi_, fw_: array [0 .. 1000]
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 // зачем тут begin и соответственно внизу предпоследний end
clrscr; // сперва выводим на экран сообщения, и не давая прочитать пользователю, сразу стираем ???? - бред.
assign(ff, 'n.txt');
rewrite(ff);
tmax := 100.0;
dt := 0.001;
for i := 1 to n do
begin
fi_[0] := 0.0; // 100 раз присвоить fi_[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); //в цикл, если надо, не попало, + см выше по ф-ции difurI
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.
Сейчас нет времени разбираться с самим алгоритмом, но код ужасный и бессмысленный. Поэтому надо не искать
, а разобраться нормально с кодом.
PS скачай Notepad++ и выбери Синтаксис-Pascal. Он тебе покажет, что к чему относиться.
PSPS Кстати, забыл самое главное:
Код:
tmax := 100.0;
dt := 0.001;
...
t := 0.0;
...
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;
while (t <= tmax) do - будет бесконечным!!!
Т.к. t наращивается вне цикла!!!
dt := 0.001;
...
t := 0.0;
...
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;
while (t <= tmax) do - будет бесконечным!!!
Т.к. t наращивается вне цикла!!!
Надоел мне Java c MySQL, решил я старый добрый паскаль вспомнить. Тока я делфи 2010 открыл.
Так вот, не вдаваясь в подробности, берем первый цикл:
при начальных параметрах, указанных автором, получаем.
При n = 30:
k1 = -3,0880777845e+85
k2 = -3,0896218234e+85
k3 = -3,0896225954e+85
k41 = -2,6670760186e+50
fi_[30] = -3,0896223382e+85
fw_[30] = -2,6657429248e+50
n =31 : k1 = -8,2361388884e+134
n =32 : k1 =-7,8963262308e+214
ну а дальше переполнение... (Real 5*10^-324 .. 1.7*10^308).
если взять все типы Extended (3.6*10^-4951 .. 1.1*10^4932-1).
То это поможет ещё где-то на 10 циклов и опять переполнение.
Так что надо разбираться в логике: Методе Рунге-Кутты (Законе Киргофа ... Коши)
Так вот, не вдаваясь в подробности, берем первый цикл:
Код:
while (n < 100) do
begin
// почему до 100, а не до 1000, если массивы объявлены так: fi_, fw_: array [0 .. 1000]
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;
begin
// почему до 100, а не до 1000, если массивы объявлены так: fi_, fw_: array [0 .. 1000]
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;
При n = 30:
k1 = -3,0880777845e+85
k2 = -3,0896218234e+85
k3 = -3,0896225954e+85
k41 = -2,6670760186e+50
fi_[30] = -3,0896223382e+85
fw_[30] = -2,6657429248e+50
n =31 : k1 = -8,2361388884e+134
n =32 : k1 =-7,8963262308e+214
ну а дальше переполнение... (Real 5*10^-324 .. 1.7*10^308).
если взять все типы Extended (3.6*10^-4951 .. 1.1*10^4932-1).
То это поможет ещё где-то на 10 циклов и опять переполнение.
Так что надо разбираться в логике: Методе Рунге-Кутты (Законе Киргофа ... Коши)