Метод Рунге-Кутты 4го порядка для уравнений двигателя постоянного тока по второму закону Киргофа

@Станиславиус · Регистрация: 23.03.2013

Студворк — интернет-сервис помощи студентам

Здравствуйте. Нужен взгляд человека понимающего 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):
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))

Собственно осуществить метод вышло на Python'e.

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
# -*- coding: UTF-8 -*-
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams['font.family'] = 'fantasy'
mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial'
from math import *
 
# Решение уравнений методом РК4
 
 
starttime = 0                                                                    #Начальное время
endtime = 10                                                                    #Конечное время
dtime = 0.1
 
# Определение количества циклов
 
ct=0
ctmax = (endtime-starttime)/dtime + 1
 
#Данные двигателя
 
U = 70.
R = 0.031
Cef = 1.2
Cmf = 1.2
L = 0.01
J = 13.
Mc = 148.
I = [0]
w1 = [0]
time1 = [0]
 
while starttime <= endtime:
    print ct
 
 
    #F1 = (U-I*R-w1*Cef)/L
    #F2 = (I*Cmf-Mc)/J
 
 
 
    k11 = dtime*((U-I[ct]*R-w1[ct]*Cef)/L)
    k21 = dtime*((U-I[ct]*R-w1[ct]*Cef)/L) + k11/2.0
    k31 = dtime*((U-I[ct]*R-w1[ct]*Cef)/L) + k21/2.0
    k41 = dtime*((U-I[ct]*R-w1[ct]*Cef)/L) + k31/2.0
 
    k12 = dtime*((I[ct]*Cmf-Mc)/J)
    k22 = dtime*((I[ct]*Cmf-Mc)/J) + k12/2.0
    k32 = dtime*((I[ct]*Cmf-Mc)/J) + k22/2.0
    k42 = dtime*((I[ct]*Cmf-Mc)/J) + k32/2.0
 
    I.append(I[ct]+1/6.0*(k11+k21*2+k31*2+k41))
    w1.append(w1[ct]+1/6.0*(k12+k22*2+k32*2+k42))
    time1.append(starttime)
 
    print ct
    ct=ct+1
    starttime = starttime+dtime
 
plt.figure(1)
#plt.subplot(211)
plt.title(u'График тока')
plt.plot(time1, I)
 
 
plt.figure(2)
#plt.subplot(212)
plt.title(u'График скорости')
plt.plot(time1, w1)
plt.show()

Теперь пытаюсь осуществить этот метод на Pascal. Выходит мягко говоря не очень, поэтому и прошу вашей помощи.
Мои потуги:

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
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(fi_,fw_,Ua,Ra,La,Cef:real):real;
 
begin
 
difurI:=(Ua-fi_*Ra-Cef*fw_)/La;
 
end;
 
function difurW(fi_,Js,Mc,Cmf: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<999) 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
 
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.

Новые блоги и статьи Все статьи Все блоги /
Как я обхитрил таблицу Word Alexander-7 21.03.2026 Когда мигает курсор у внешнего края таблицы, и нам надо перейти на новую строку, а при нажатии Enter создается новый ряд таблицы с ячейками, то мы вместо нервных нажатий Энтеров мы пишем любые буквы. . .	Krabik - рыболовный бот для WoW 3.3.5a AmbA 21.03.2026 без регистрации и смс. Это не торговля, приложение не содержит рекламы. Выполняет свою непосредственную задачу - автоматизацию рыбалки в WoW - и ничего более. Однако если админы будут против -. . .	Программный отбор значений справочника Maks 21.03.2026 Установка программного отбора значений справочника "Сотрудники" из модуля формы документа. В качестве фильтра для отбора служит предопределенное значение перечислений. Процедура. . .	Переходник USB-CAN-GPIO Eddy_Em 20.03.2026 Достаточно давно на работе возникла необходимость в переходнике CAN-USB с гальваноразвязкой, оный и был разработан. Однако, все меня терзала совесть, что аж 48-ногий МК используется так тупо: просто. . .
Оттенки серого Argus19 18.03.2026 Оттенки серого Нашёл в интернете 3 прекрасных модуля: Модуль класса открытия диалога открытия/ сохранения файла на Win32 API; Модуль класса быстрого перекодирования цветного изображения в оттенки. . .	SDL3 для Desktop (MinGW): Рисуем цветные прямоугольники с помощью рисовальщика SDL3 на Си и C++ 8Observer8 17.03.2026 Содержание блога Финальные проекты на Си и на C++: finish-rectangles-sdl3-c. zip finish-rectangles-sdl3-cpp. zip	Символические и жёсткие ссылки в Linux. algri14 15.03.2026 Существует два типа ссылок — символические и жёсткие. Ссылка в Linux — это запись в каталоге, которая может указывать либо на inode «файла-ИСТОЧНИКА», тогда это будет «жёсткая ссылка» (hard link),. . .	[Owen Logic] Поддержание уровня воды в резервуаре количеством включённых насосов: моделирование и выбор регулятора ФедосеевПавел 14.03.2026 Поддержание уровня воды в резервуаре количеством включённых насосов: моделирование и выбор регулятора ВВЕДЕНИЕ Выполняя задание на управление насосной группой заполнения резервуара,. . .