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

@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 **[ТС]**

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

Новые блоги и статьи Все статьи Все блоги /
Кто-нибудь знает, где можно бесплатно получить настольный компьютер или ноутбук? США. Programma_Boinc 26.12.2025 Нашел на реддите интересную статью под названием Anyone know where to get a free Desktop or Laptop? Ниже её машинный перевод. После долгих разбирательств я наконец-то вернула себе. . .	Thinkpad X220 Tablet — это лучший бюджетный ноутбук для учёбы, точка. Programma_Boinc 23.12.2025 Рецензия / Мнение/ Перевод Нашел на реддите интересную статью под названием The Thinkpad X220 Tablet is the best budget school laptop period . Ниже её машинный перевод. Thinkpad X220 Tablet —. . .	PhpStorm 2025.3: WSL Terminal всегда стартует в ~ and_y87 14.12.2025 PhpStorm 2025. 3: WSL Terminal всегда стартует в ~ (home), игнорируя директорию проекта Симптом: После обновления до PhpStorm 2025. 3 встроенный терминал WSL открывается в домашней директории. . .	Как объединить две одинаковые БД Access с разными данными VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.	Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro
Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .	Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .

@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