Как решить систему ДУ для системы маятников методом Рунге-Кутта?

@patzo · Регистрация: 04.12.2012

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

Помогите, пожалуйста!
Не работает программа. Суть задания - решить систему ДУ для системы маятников. Не понимаю, почему не работает - насколько сумел разобраться - что-то не так с заданием функции, но что конкретно - понять не могу(

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
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
""" Cистема маятников """
from numpy import * 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import time as tm
from math import *
from scipy import linspace, array
 
starttime = 0                                                                    #Начальное время
endtime = 10                                                                    #Конечное время
tmax = 100 #задание количества элементарных отрезков
# Определение количества циклов 
t = starttime
dt = (endtime-starttime)/tmax + 1
#Данные системы маятников
m2 = 2 # масса нижнего маятника
m1 = 2*m2 # масса верхнего маятника
g = 9.8 # ускорение свободного падения
l = 5 # длина нитей
alfa1 = [] # частота колебаний верхнего маятника 
alfa2 = [] # частота колебаний нижнего маятника
time = [0] # время отсчёта
p1 = [] # импульс верхнего маятника
p2 = [] # импульс нижнего маятника
alfa1[0] = 0
alfa2[0] = 0
p1[0] = 0
p2[0] = 0
 
 
 
 
#def F(Z, t)
def F1(alfa1, alfa2, p1, p2, t):
    
    F1 = (p1-p2*cos(alfa1-alfa2))/(m1*l*(1+m2/m1*sin(alfa1-alfa2)**2))   #Функция угла 1-го маятника
    
    return [alfa1, alfa2, p1, p2, t]
    
def F2(alfa1, alfa2, p1, p2, t):  
    
    F2 = (p2*(1+m2/m1)-p1*m2/m1*cos(alfa1-alfa2))/(m1*l**2*(1+m2/m1*sin(alfa1-alfa2)**2)) #Функция угла 2-го маятника
    
    return [alfa1, alfa2, p1, p2, t]
    
def F3(alfa1, alfa2, p1, p2, t):
    
    F3 = -m1*(1+m2/m1)*g*l*sin(alfa1) - p1*p2*sin(alfa1-alfa2)/(l**2*(m1+m2*sin(alfa1-alfa2)**2)) + (p1**2*m2*l**2-2*p1*p2*l**2*cos(alfa1-alfa2)+p2**2*(m1+m2)*l**2)*sin(2*(alfa1-alfa2))/(2*l**4*(m1+m2*sin(alfa1-alfa2)**2)**2) #Функция импульса 1-го маятника
    
    return [alfa1, alfa2, p1, p2, t]
    
def F4(alfa1, alfa2, p1, p2, t):
    
    F4 = -m1*m2/m1*g*l*sin(alfa2) + p1*p2*sin(alfa1-alfa2)/(l**2*(m1+m2*sin(alfa1-alfa2)**2)) - (p1**2*m2*l**2-2*p1*p2*l**2*cos(alfa1-alfa2)+p2**2*(m1+m2)*l**2)*sin(2*(alfa1-alfa2))/(2*l**4*(m1+m2*sin(alfa1-alfa2)**2)**2)  #Функция импульса 2-го маятника
 
    return [alfa1, alfa2, p1, p2, t]
 
#while starttime <= endtime
 
for i in range(0, tmax-1):
     #  коэффициенты для функции alfa1(t), alfa2(t), p1(t), p2(t)
    K1 = F1(t, p1[i], p2[i], alfa1[i], alfa2[i])  #- это для уравнения alfa1'(t)
    L1 = F2(t, p1[i], p2[i], alfa1[i], alfa2[i])  #- это для уравнения alfa2'(t)
    M1 = F3(t, p1[i], p2[i], alfa1[i], alfa2[i])  #- это для уравнения p1'(t)
    N1 = F4(t, p1[i], p2[i], alfa1[i], alfa2[i])  #- это для уравнения p2'(t)
 
    K2 = F1(t +dt/2.0, p1[i] + M1/2.0, p2[i] + N1/2.0, alfa1[i] + K1/2.0, alfa2[i] + L1/2.0)
    L2 = F2(t +dt/2.0, p1[i] + M1/2.0, p2[i] + N1/2.0, alfa1[i] + K1/2.0, alfa2[i] + L1/2.0)
    M2 = F3(t +dt/2.0, p1[i] + M1/2.0, p2[i] + N1/2.0, alfa1[i] + K1/2.0, alfa2[i] + L1/2.0)
    N2 = F4(t +dt/2.0, p1[i] + M1/2.0, p2[i] + N1/2.0, alfa1[i] + K1/2.0, alfa2[i] + L1/2.0)
    
    K3 = F1(t +dt/2.0, p1[i] + M2/2.0, p2[i] + N2/2.0, alfa1[i] + K2/2.0, alfa2[i] + L2/2.0)
    L3 = F2(t +dt/2.0, p1[i] + M2/2.0, p2[i] + N2/2.0, alfa1[i] + K2/2.0, alfa2[i] + L2/2.0)
    M3 = F3(t +dt/2.0, p1[i] + M2/2.0, p2[i] + N2/2.0, alfa1[i] + K2/2.0, alfa2[i]+ L2/2.0)
    N3 = F4(t +dt/2.0, p1[i] + M2/2.0, p2[i] + N2/2.0, alfa1[i] + K2/2.0, alfa2[i] + L2/2.0)
    
    K4 = F1(t +dt/2.0, p1[i] + M3/2.0, p2[i] + N3/2.0, alfa1[i] + K3/2.0, alfa2[i] + L3/2.0)
    L4 = F2(t +dt/2.0, p1[i] + M3/2.0, p2[i] + N3/2.0, alfa1[i] + K3/2.0, alfa2[i] + L3/2.0)
    M4 = F3(t +dt/2.0, p1[i] + M3/2.0, p2[i] + N3/2.0, alfa1[i] + K3/2.0, alfa2[i] + L3/2.0)
    N4 = F4(t +dt/2.0, p1[i] + M3/2.0, p2[i] + N3/2.0, alfa1[i] + K3/2.0, alfa2[i] + L3/2.0)
    
    #Вычисляем значения функций в i+1 узле:
        
    alfa1[i+1] = alfa1[i] + dt/6*(K1 + 2*K2 + 2*K3 + K4)
    alfa2[i+1] = alfa2[i] + dt/6*(L1 + 2*L2 + 2*L3 + L4)
    p1[i+1] = p1[i] + dt/6*(M1 + 2*M2 + 2*M3 + M4)
    p2[i+1] = p2[i] + dt/6*(N1 + 2*N2 + 2*N3 + N4)
 
    t = t + dt
    
print ('alfa1=', alfa1[i+1])
print ('alfa2=', alfa2[i+1])
print ('p1=', p1[i+1])
print ('p2=', p2[i+1])

Часть библиотек - лишние - копировал с большой программы...

@Wolkodav · 24.10.2013, 22:21

patzo, что пишет?

@patzo · 24.10.2013, 23:48 **[ТС]**

Python
1
2
3
4
5
6
7
8
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Python27\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 487, in runfile
    execfile(filename, namespace)
  File "C:\Users\patzo\Documents\Python\Hometask 1.py", line 248, in <module>
    K2 = F1(t +dt/2.0, p1[i] + M1/2.0, p2[i] + N1/2.0, alfa1[i] + K1/2.0, alfa2[i] + L1/2.0)
TypeError: unsupported operand type(s) for /: 'list' and 'float'
>>>

Если же написать везде

Python
1
return array([alfa1, ... ])

то ругается на то, что мол должен быть один аргумент.

@dondublon · 28.10.2013, 14:43

Не по теме:

Он метод Рунге-Кутты

@patzo · 28.10.2013, 21:17 **[ТС]**

Не по теме:
Его и так и так называют.

Никто не поможет?(

@Wolkodav · 28.10.2013, 22:29

Не по теме:

patzo, её, уж если на то пошло.

Завтра на работе свободная минуточка будет, гляну.

accept · 28.10.2013, 23:05

Не по теме:

Сообщение от Wolkodav

patzo, её, уж если на то пошло.

кого её ? Кутта - это мужик

@patzo · 01.11.2013, 16:47 **[ТС]**

Никто не знает?

Новые блоги и статьи Все статьи Все блоги /
Модульный подход на примере F# DevAlt 06.03.2026 В блоге дяди Боба наткнулся на такое определение: В этой книге («Подход, основанный на вариантах использования») Ивар утверждает, что архитектура программного обеспечения — это структуры,. . .	Управление камерой с помощью скрипта OrbitControls.js на Three.js: Вращение, зум и панорамирование 8Observer8 05.03.2026 Содержание блога Финальная демка в браузере работает на Desktop и мобильных браузерах. Итоговый код: orbit-controls-threejs-js. zip. Сканируйте QR-код на мобильном. Вращайте камеру одним пальцем,. . .	SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .	SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip Сканируйте QR-код на мобильном и вы увидите, что появится джойстик для управления главным героем. . . .
Реалии Hrethgir 01.03.2026 Нет, я не закончил до сих пор симулятор. Эта задача сложнее. Не получилось уйти в плавсостав, но оно и к лучшему, возможно. Точнее получалось - но сварщиком в палубную команду, а это значит, в моём. . .	Ритм жизни kumehtar 27.02.2026 Иногда приходится жить в ритме, где дел становится всё больше, а вовлечения в происходящее — всё меньше. Плотный график не даёт вниманию закрепиться ни на одном событии. Утро начинается с быстрых,. . .	SDL3 для Web (WebAssembly): Сборка библиотек: SDL3, Box2D, FreeType, SDL3_ttf, SDL3_mixer и SDL3_image из исходников с помощью CMake и Emscripten 8Observer8 27.02.2026 Недавно вышла версия 3. 4. 2 библиотеки SDL3. На странице официальной релиза доступны исходники, готовые DLL (для x86, x64, arm64), а также библиотеки для разработки под Android, MinGW и Visual Studio. . . .	SDL3 для Web (WebAssembly): Реализация движения на Box2D v3 - трение и коллизии с повёрнутыми стенами 8Observer8 20.02.2026 Содержание блога Box2D позволяет легко создать главного героя, который не проходит сквозь стены и перемещается с заданным трением о препятствия, которые можно располагать под углом, как верхнее. . .

@Wolkodav 842 / 480 / 58 Регистрация: 18.09.2012 Сообщений: 1,688
	24.10.2013, 22:21
	patzo, что пишет? 0

@dondublon
	28.10.2013, 14:43
	Не по теме: Он метод Рунге-Кутты 0

@patzo 0 / 0 / 0 Регистрация: 04.12.2012 Сообщений: 7
	28.10.2013, 21:17 [ТС]
	Не по теме: Его и так и так называют. Никто не поможет?( 0

@Wolkodav 842 / 480 / 58 Регистрация: 18.09.2012 Сообщений: 1,688
	28.10.2013, 22:29
	Не по теме: patzo, её, уж если на то пошло. Завтра на работе свободная минуточка будет, гляну. 0

@patzo 0 / 0 / 0 Регистрация: 04.12.2012 Сообщений: 7
	01.11.2013, 16:47 [ТС]
	Никто не знает? 0