Зациклить решение системы дифф-уравнений пока есть возможные варианты значений переменных

@virus77724 · Регистрация: 22.12.2015

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

Есть программа которая решает систему дифференциальных уравнений методом Адамса Башфорта.

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
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
import scipy.integrate
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg 
from matplotlib.figure import Figure 
 
def GetF(k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16, k17, k18, k19, k20, k21, k22, k23, k24, k25, k26, k27, k28, k29, k30, k31, k32, k33, k34, k35, k36, k37, k38, k39, k40, k41, k42, k43, k44, k45, k46, k47, k48, k49, k50, k51, k52, k53, k54, k55, k56, k57, k58, k59, k60, k61, k62, k63, k64, k65, k66, k67):
    f1 = lambda t, yList : 0.011*k1*yList[10] + 0.019*k10*yList[10] + 0.05*k13*yList[11] + 0.07*k22*yList[11] + 0.04*k23*yList[9] + 0.04*k31*yList[9] + 0.11*k32*yList[8] + 0.05*k39*yList[8] + 0.07*k40*yList[7] + 0.22*k46*yList[7] + 0.47*k47*yList[6] + 0.03*k52*yList[6] + 0.61*k53*yList[5] + 0.08*k57*yList[5] + 0.3*k58*yList[4] + 0.16*k61*yList[4] + 0.8*k62*yList[3] + 0.11*k64*yList[3] + 0.37*k65*yList[2] + 0.63*k66*yList[3] + 1*k67*yList[1]
    f2 = lambda t, yList : 0.15*k2*yList[10] + 0.01*k9*yList[10] + 0.09*k14*yList[11] + 0.07*k21*yList[11] + 0.07*k24*yList[9] + 0.08*k30*yList[9] + 0.09*k33*yList[8] + 0.12*k38*yList[8] + 0.11*k41*yList[7] + 0.06*k45*yList[7] + 0.23*k48*yList[6] + 0.05*k51*yList[6] + 0.18*k54*yList[5] + 0.12*k56*yList[5] + 0.35*k59*yList[4] + 0.19*k60*yList[4] + 0.09*k63*yList[3] + 0.09*k63*yList[3] + 0.63*k66*yList[2] + 0.37*k65*yList[3] - 1*k67*yList[1]
    f3 = lambda t, yList : 0.1*k3*yList[10] + 0.03*k8*yList[10] + 0.12*k15*yList[11] + 0.09*k20*yList[11] + 0.17*k25*yList[9] + 0.08*k29*yList[9] + 0.33*k34*yList[8] + 0.05*k37*yList[8] + 0.27*k42*yList[7] + 0.18*k44*yList[7] + 0.14*k49*yList[6] + 0.08*k50*yList[6] + 0.01*k55*yList[5] + 0.01*k55*yList[5] + 0.19*k60*yList[4] + 0.35*k59*yList[4] + 0.11*k64*yList[3] + 0.8*k62*yList[3] - 0.37*k65*yList[3] - 0.63*k66*yList[3]
    f4 = lambda t, yList : 0.07*k4*yList[10] + 0.04*k7*yList[10] + 0.07*k16*yList[11] + 0.11*k19*yList[11] + 0.22*k26*yList[9] + 0.12*k28*yList[9] + 0.17*k35*yList[8] + 0.08*k36*yList[8] + 0.09*k43*yList[7] + 0.09*k43*yList[7] + 0.08*k50*yList[6] + 0.14*k49*yList[6] + 0.12*k56*yList[5] + 0.18*k54*yList[5] + 0.16*k61*yList[4] + 0.3*k58*yList[4] - 0.11*k64*yList[3] - 0.09*k63*yList[3] - 0.8*k62*yList[3]
    f5 = lambda t, yList : 0.31*k5*yList[10] + 0.06*k6*yList[10] + 0.2*k17*yList[11] + 0.13*k18*yList[11] + 0.18*k27*yList[9] + 0.18*k27*yList[9] + 0.08*k36*yList[8] + 0.17*k35*yList[8] + 0.18*k44*yList[7] + 0.27*k42*yList[7] + 0.05*k51*yList[6] + 0.23*k48*yList[6] + 0.08*k57*yList[5] + 0.61*k53*yList[5] - 0.3*k58*yList[4] - 0.35*k59*yList[4] - 0.19*k60*yList[4] - 0.16*k61*yList[4]
    f6 = lambda t, yList : 0.06*k6*yList[10] + 0.31*k5*yList[10] + 0.13*k18*yList[11] + 0.2*k17*yList[11] + 0.12*k28*yList[9] + 0.22*k26*yList[9] + 0.05*k37*yList[8] + 0.33*k34*yList[8] + 0.06*k45*yList[7] + 0.11*k41*yList[7] + 0.03*k52*yList[6] + 0.47*k47*yList[6] - 0.61*k53*yList[5] - 0.18*k54*yList[5] - 0.01*k55*yList[5] - 0.12*k56*yList[5] - 0.08*k57*yList[5]
    f7 = lambda t, yList : 0.04*k7*yList[10] + 0.07*k4*yList[10] + 0.11*k19*yList[11] + 0.07*k16*yList[11] + 0.08*k29*yList[9] + 0.17*k25*yList[9] + 0.12*k38*yList[8] + 0.09*k33*yList[8] + 0.22*k46*yList[7] + 0.07*k40*yList[7] - 0.47*k47*yList[6] - 0.23*k48*yList[6] - 0.14*k49*yList[6] - 0.08*k50*yList[6] - 0.05*k51*yList[6] - 0.03*k52*yList[6]
    f8 = lambda t, yList : 0.03*k8*yList[10] + 0.1*k3*yList[10] + 0.09*k20*yList[11] + 0.12*k15*yList[11] + 0.08*k30*yList[9] + 0.07*k24*yList[9] + 0.05*k39*yList[8] + 0.11*k32*yList[8] - 0.07*k40*yList[7] - 0.11*k41*yList[7] - 0.27*k42*yList[7] - 0.09*k43*yList[7] - 0.18*k44*yList[7] - 0.06*k45*yList[7] - 0.22*k46*yList[7]
    f9 = lambda t, yList : 0.02*k9*yList[10] + 0.15*k2*yList[10] + 0.07*k21*yList[11] + 0.09*k14*yList[11] + 0.04*k31*yList[9] + 0.04*k23*yList[9] - 0.11*k32*yList[8] - 0.09*k33*yList[8] - 0.33*k34*yList[8] - 0.17*k35*yList[8] - 0.08*k36*yList[8] - 0.05*k37*yList[8] - 0.12*k38*yList[8] - 0.05*k39*yList[8]
    f10 = lambda t, yList : 0.01*k9*yList[10] + 0.15*k2*yList[10] + 0.07*k21*yList[11] + 0.09*k14*yList[11] - 0.04*k23*yList[9] - 0.07*k24*yList[9] - 0.17*k25*yList[9] - 0.22*k26*yList[9] - 0.18*k27*yList[9] - 0.12*k28*yList[9] - 0.08*k29*yList[9] - 0.08*k30*yList[9] - 0.04*k31*yList[9]
    f11 = lambda t, yList : 0.12*k12*yList[10] - 0.011*k1*yList[10] - 0.15*k2*yList[10] - 0.1*k3*yList[10] - 0.07*k4*yList[10] - 0.31*k5*yList[10] - 0.06*k6*yList[10] - 0.04*k7*yList[10] - 0.03*k8*yList[10] - 0.02*k9*yList[10] - 0.019*k10*yList[10]
    f12 = lambda t, yList : 0.18*k11*yList[11] - 0.05*k13*yList[11] - 0.09*k14*yList[11] - 0.12*k15*yList[11] - 0.07*k16*yList[11] - 0.2*k17*yList[11] - 0.13*k18*yList[11] - 0.11*k19*yList[11] - 0.09*k20*yList[11] - 0.07*k21*yList[11] - 0.07*k22*yList[11]
    
    return lambda t, yList : np.asarray([f1(t, yList), f2(t, yList), f3(t, yList), f4(t, yList), f5(t, yList), f6(t, yList), f7(t, yList), f8(t, yList), f9(t, yList), f10(t, yList), f11(t, yList), f12(t, yList)])
 
def calc_k():
    t0 = float(t0_entry.get())
    tLast = float(tLast_entry.get())
    h = float(h_entry.get())
    y10 = float(y10_entry.get())
    y0=y1=y2=y3=y4=y5=y6=y7=y8=y9=y11=0
    
    k1 = 1.6; k2 = 0.5; k3 = 2.3; k4 = 1.7; k5 = 3.8; k6 = 5.9; k7 = 1.1; k8 = 2.2; k9 = 0.4; k10 = 3.1; k11 = 6.8; k12 = 12.8
    k13 = 2.3; k14 = 1.4; k15 = 1.5; k16 = 0.6; k17 = 6.8; k18 = 2.9; k19 = 1.2; k20 = 0.3; k21 = 3.7; k22 = 2.1; k23 = 0.4; k24 = 1.5
    k25 = 1.3; k26 = 0.2; k27 = 0.1; k28 = 5.8; k29 = 8.7; k30 = 2.9; k31 = 1.4; k32 = 3.6; k33 = 1.4; k34 = 1.5; k35 = 5.9; k36 = 9.8
    k37 = 3.6; k38 = 1.3; k39 = 0.1; k40 = 0.4; k41 = 2.2; k42 = 4.3; k43 = 6.8; k44 = 1.3; k45 = 1.5; k46 = 0.2; k47 = 3.8; k48 = 4.4
    k49 = 0.3; k50 = 1.8; k51 = 2.2; k52 = 2.1; k53 = 0.19; k54 = 3.4; k55 = 9.6; k56 = 1.2; k57 = 0.7; k58 = 1.8; k59 = 0.6
    k60 = 1.4; k61 = 2.1; k62 = 0.2; k63 = 0.8; k64 = 0.3; k65 = 0.3; k66 = 0.2; k67 = 0.3
 
# Функция правых частей в векторной форме
 
    f = GetF(k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16, k17, k18, k19, k20, k21, k22, k23, k24, k25, k26, k27, k28, k29, k30, k31, k32, k33, k34, k35, k36, k37, k38, k39, k40, k41, k42, k43, k44, k45, k46, k47, k48, k49, k50, k51, k52, k53, k54, k55, k56, k57, k58, k59, k60, k61, k62, k63, k64, k65, k66, k67)
    rk = scipy.integrate.RK45(
        fun = f,
        t0 = t0, 
        y0 = [y0, y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11],
        t_bound = tLast)
 
    rkPoints = []
    time = []
 
# Метод Рунге-Кутта
    
    for i in range(4):
        try:
            rkPoints.append(np.asarray(rk.y))
            time.append(rk.t)
            rk.step()
        except Exception as ex:
            break
 
# Метод Адамса-Башфорта
    dt = h
    for i in range( int((tLast-t0)/h) ):
        newPoint = rkPoints[-1] + h/24 * (55*f(time[-1], rkPoints[-1]) -
             59*f(time[-2], rkPoints[-2]) +
             37*f(time[-3], rkPoints[-3]) -
              9*f(time[-4], rkPoints[-4]) )
        rkPoints.append(newPoint)
        time.append(time[-1] + dt)
    
# Построение графика    
 
    fig = plt.figure(figsize=(16, 9), dpi=200)
    a = fig.add_subplot(111)
    a.set_xlabel('Время, с')
    a.set_ylabel('Концентрация, Моль')
    plt.plot(time, rkPoints)
    a.legend(("Метан", "Этан","Пропан","Бутан","Пентан", "Гексан","Гептан","Октан","Нонан", "Декан","Ундекан","Изо-Ундекан"))
    
# Сохранение результатов отдельно в корне программы
 
    plt.savefig("Graph.png")
    f = open('Rezultat.txt', 'w')
    print (rkPoints[i], file=f) 
    f.close()    
        
    canvas = FigureCanvasTkAgg(fig) 
    canvas.get_tk_widget().pack() 
    canvas.draw() 
              
# создание экземпляра класса Tk, отвечающего за создание окон
root = tkinter.Tk ()
# определение заголовка окна
root. title ("Решения системы дифференциальных уравнений методами Рунге-Кутта и Адамса Башфорта")
frame = tkinter.Frame(root)
frame. pack ()
#########################################################################################################################
# создание окна ввода начальной концентрация для y11
y10_entry = tkinter.Entry (frame, width=10)
y10_entry.insert(0, "0.3")
y10_entry. grid (row=1, column=9,pady=5)
y10_label = tkinter.Label (frame, text="Начальная концентрация Ундекана: ")
y10_label. grid (row=1, column=8,pady=5)
######################################################################################################################
# создание окна ввода конечного момента времени
tLast_entry = tkinter.Entry (frame, width=10)
tLast_entry.insert(0, "1")
tLast_entry. grid (row=2, column=9,pady=5)
tLast_label = tkinter.Label (frame, text="Конечный момент времени (правая граница): ")
tLast_label. grid (row=2, column=8,pady=5)
# создание окна ввода значения парамметра t при котором получены y1, y2, y3, y4:
t0_entry = tkinter.Entry (frame, width=10)
t0_entry.insert(0, "0")
t0_entry. grid (row=3, column=9,pady=5)
t0_lebel = tkinter.Label (frame, text="Значение парамметра t  ")
t0_lebel. grid (row=3, column=8,pady=5)
# создание окна ввода величины точности интегрирования)
h_entry = tkinter.Entry (frame, width=10)
h_entry.insert(0, "0.001")
h_entry. grid (row=4, column=9,pady=5)
h_label = tkinter.Label (frame, text="Размер разбиения: ")
h_label. grid (row=4, column=8,pady=5)
##########################################################################################################################
# создание кнопки решения
eval_button = tkinter.Button (frame,bg='coral', text="Решить", width=10,command=calc_k)
eval_button. grid (row=38, column=8,padx=5,pady=5)
# создание кнопки закрытия приложения
exit_button = tkinter.Button (frame, bg='coral', text="Выход", width=10,command=root. destroy)
exit_button. grid (row=38, column=9,padx=5,pady=5)
# создание окна
root. mainloop ()

Нужно переменным "k" дать не строгое значение, а диапазон от 0.1 до 1.0. Соответственно создать цикл чтобы решать систему со всеми возможными комбинациями возможных "k".
(67*10)*(67*10)=448900 возможных комбинаций "k", а значит и повторов выполнения

@dondublon · 24.04.2019, 14:10

virus77724, ужас какой.
Откройте для себя массивы, раз уж вы залезли в эти дебри.

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

Нужно переменным "k" дать не строгое значение, а диапазон от 0.1 до 1.0. Соответственно создать цикл чтобы решать систему со всеми возможными комбинациями возможных "k".

А в чём проблема?

Новые блоги и статьи Все статьи Все блоги /
AkelPad-скрипты, структуры, и немного лирики.. testuser2 05.04.2026 Такая программа, как AkelPad существует уже давно, и также давно существуют скрипты под нее. Тем не менее, прога живет, периодически что-то не спеша дополняется, улучшается. Что меня в первую очередь. . .	Отображение реквизитов в документе по условию и контроль их заполнения Maks 04.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "ПланированиеСпецтехники", разработанного в конфигурации КА2. Данный документ берёт данные из другого нетипового документа. . .	Фото всей Земли с борта корабля Orion миссии Artemis II kumehtar 04.04.2026 Это первое подобное фото сделанное человеком за 50 лет. Снимок называют новым вариантом легендарной фотографии «The Blue Marble» 1972 года, сделанной с борта корабля «Аполлон-17». Новое фото. . .	Вывод диалогового окна перед закрытием, если документ не проведён Maks 04.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "СписаниеМатериалов", разработанного в конфигурации КА2. Задача: реализовать программный контроль на предмет проведения документа. . .
Программный контроль заполнения реквизитов табличной части документа Maks 02.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "СписаниеМатериалов", разработанного в конфигурации КА2. Задача: 1. Реализовать контроль заполнения реквизита. . .	wmic не является внутренней или внешней командой Maks 02.04.2026 Решение: DISM / Online / Add-Capability / CapabilityName:WMIC~~~~ Отсюда: https:/ / winitpro. ru/ index. php/ 2025/ 02/ 14/ komanda-wmic-ne-naydena/	Программная установка даты и запрет ее изменения Maks 02.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "СписаниеМатериалов", разработанного в конфигурации КА2. Задача: при создании документов установить период списания автоматически. . .	Вывод данных в справочнике через динамический список Maks 01.04.2026 Реализация из решения ниже выполнена на примере нетипового справочника "Спецтехника" разработанного в конфигурации КА2. Задача: вывести данные из ТЧ нетипового документа. . .