{"id":185,"date":"2022-01-11T19:30:06","date_gmt":"2022-01-11T18:30:06","guid":{"rendered":"http:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/?p=185"},"modified":"2022-02-06T15:39:13","modified_gmt":"2022-02-06T14:39:13","slug":"simulace-ode-v-pythonu","status":"publish","type":"post","link":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/2022\/01\/11\/simulace-ode-v-pythonu\/","title":{"rendered":"Simulace ODE v Pythonu"},"content":{"rendered":"\n

Dnes se pod\u00edv\u00e1me na problematiku simulace oby\u010dejn\u00e9 diferenci\u00e1ln\u00ed rovnice ODE v pythonu. Budeme vyu\u017e\u00edvat modul scipy, kter\u00fd obsahuje mnoho algortim\u016f pro numerickou matematiku. Pro \u0159e\u0161en\u00ed dynamick\u00e9ho syst\u00e9m\u016f (ten je pops\u00e1n diferenci\u00e1ln\u00edmi rovnicemi) pou\u017eijeme numerickou metodu. Rozd\u00edl mezi numerick\u00fdm a analytick\u00e9m \u0159e\u0161en\u00edm je ten, \u017ee analytick\u00e9 \u0159e\u0161en\u00ed je zcela p\u0159esn\u00e9, ale lze takto vy\u0159e\u0161it jen omezenou sadu \u00faloh. Naopak numerika je velmi pou\u017e\u00edvan\u00e1 v in\u017een\u00fdrsk\u00e9 praxi. <\/p>\n\n\n\n

Jak vypad\u00e1 takov\u00e1 ODE?<\/p>\n\n\n\n

\\large\\frac{dx}{dt}= ax <\/pre><\/div>\n\n\n\n
Tak nap\u0159\u00edklad v\u00fd\u0161e m\u00e1me ODE prvn\u00edho \u0159\u00e1du. Na lev\u00e9 stran\u011b m\u00e1me derivaci x<\/em> podle \u010dasu.  Pro tuto jednoduchou ODE je zn\u00e1m\u00e9 analytick\u00e9 \u0159e\u0161en\u00ed. \u0158e\u0161en\u00edm ODE je funkce:<\/p>\n\n\n\n
\\large x(t)= e^{at}x(0)<\/pre><\/div>\n\n\n\n
Take se m\u016f\u017eeme setkat s term\u00ednem „\u0159e\u0161en\u00ed po\u010d\u00e1te\u010dn\u00edho probl\u00e9mu“. M\u016fs\u00edme tedy zn\u00e1t po\u010d\u00e1te\u010dn\u00ed podm\u00ednku x(0)<\/em>. <\/p>\n\n\n\n
My budeme \u0159e\u0161it rovnice numerickou metodou. K tomu vyu\u017eijeme modul scipy a funkci solve_ipv. V prv\u00e9 \u0159ad\u011b mus\u00edme nadefinovat funkci, kterou bude numerick\u00e1 metoda volat a z\u00edsk\u00e1vat tak \u0159e\u0161en\u00ed v ka\u017ed\u00e9m v\u00fdpa\u010detn\u00edm kroku. D\u00e1le interval na jak\u00e9m budeme \u0159e\u0161en\u00ed cht\u00edt nap\u0159. \u010das simulace. Bude pot\u0159eba zadata po\u010d\u00e1te\u010dn\u00ed podm\u00ednky.  Funkce solve_ipv tak\u00e9 nab\u00edz\u00ed n\u011bkolik numerick\u00fdch metod. Pro za\u010d\u00e1tek je nejvhodn\u011bj\u0161\u00ed volit RK45. Jedn\u00e1 se o metodu s adaptivn\u00edm krokem simulace. Velikost v\u00fdpo\u010detn\u00edho kroku upravuje dle odhadu chyby \u0159e\u0161en\u00ed. Konkr\u00e9tn\u011b se jedn\u00e1 o motedu z rodiny Runge\u2013Kutta\u2013Fehlberg, co\u017e je velmi dobr\u00e1 explicitn\u00ed metoda 4-5 \u0159\u00e1du. Pro n\u011bkter\u00e9 typy \u00faloh v\u0161ak vhodn\u00e1 nen\u00ed. jedn\u00e1 se zejm\u00e9na o tzv. stiff probl\u00e9my. <\/p>\n\n\n\n
Vy\u0159e\u0161\u00edme jednoduchou ODE1 (prvn\u00edho \u0159\u00e1du)<\/p>\n\n\n\n
\"\"<\/figure><\/div>\n\n\n\n
\u03c4<\/em> = 0.1 
k<\/em> = 10
u<\/em> <0-0.1> = 0 ; <0.1-1> = 1<\/p>\n\n\n\n
Vstupem tedy bude jednotkov\u00fd skok v \u010dase 0,1 s<\/em>.<\/p>\n\n\n\n
P\u0159\u00edklad \u0159e\u0161en\u00ed vykreslen v matplotlib:<\/p>\n\n\n\n
\"\"<\/figure><\/div>\n\n\n\n
D\u00e1le vy\u0159e\u0161\u00edme rovnici, kter\u00e1 je asi nejpou\u017e\u00edvan\u011bj\u0161\u00ed ve \u0161koln\u00edm prost\u0159ed\u00ed. Jedn\u00e1 se o hmotu s pru\u017einou a tlumi\u010dem. Tento syst\u00e9m m\u00e1 i pom\u011brn\u011b dobr\u00e9 praktick\u00e9 vyu\u017eit\u00ed nap\u0159. nejednodu\u0161\u0161\u00ed model zav\u011b\u0161en\u00ed n\u00e1prav automobilu. <\/p>\n\n\n\n
\"\"<\/figure><\/div>\n\n\n\n
\\large \\ddot{x}m+\\dot{x}b+ kx = F(t) <\/pre><\/div>\n\n\n\n
Setrva\u010dn\u00e1, tlum\u00edc\u00ed a potencion\u00e1ln\u00ed (s\u00edla pru\u017einy) je v tomto p\u0159\u00edpad\u011b rovna bud\u00edc\u00ed s\u00edle. Tentokr\u00e1t jsme pou\u017eili tzv. Newton\u016fv z\u00e1pis. Jedna te\u010dka nad x<\/em> znam\u00e9n\u00e1 prvn\u00ed derivaci podle \u010dasu atd.<\/p>\n\n\n\n
P\u0159i prvn\u00edm pohledu zjist\u00edme, \u017ee se jedn\u00e1 o ODE druh\u00e9ho \u0159\u00e1du, co\u017e je typick\u00e9 pro mechanick\u00e9 syst\u00e9my. Numerick\u00e9 \u0159e\u0161i\u010de \u0159e\u0161\u00ed soustavy ODE prvn\u00edho \u0159\u00e1du a je tedy nutn\u00e9 na\u0161i rovnici upravit. Tato \u00faprava je mo\u017en\u00e1 v\u017edy a tedy rovnice vy\u0161\u0161\u00edch \u0159\u00e1d\u016f m\u016f\u017eeme p\u0159ev\u00e9st na soustavu rovnic \u0159\u00e1du prvn\u00edho.<\/p>\n\n\n\n
Upraven\u00e1 ovnice vhodn\u00e1 pro simulaci bude vypadat takto:<\/p>\n\n\n\n
\\large \\dot{x}=v \\\\ \\dot{v} = \\frac{1}{m}(F(t) - vb - xk)<\/pre><\/div>\n\n\n\n
Takto upraven\u00e1 rovnice je u\u017e vhodn\u00e1 k simulaci. Jedn\u00e1 se o soustavu ODE, kde m\u00e1me prvn\u00ed derivace izolovan\u00e9 na jedn\u00ed stran\u011b rovnic.<\/p>\n\n\n\n
\"\"<\/figure><\/div>\n\n\n\n
\"\"<\/figure><\/div>\n\n\n\n
\"\"<\/figure><\/div>\n\n\n
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
<\/div><\/td>
import<\/span> numpy as<\/span> np
\nfrom<\/span> scipy import<\/span> integrate
\nimport<\/span> matplotlib.pyplot<\/span> as<\/span> plt
\n
\n
\ndef<\/span> ode1(<\/span>t,<\/span> x)<\/span>:
\n    tau =<\/span> 0.1<\/span>
\n    k =<\/span> 10<\/span>
\n
\n    if<\/span> t <<\/span> 0.1<\/span>:
\n        u =<\/span> 0<\/span>
\n    else<\/span>:
\n        u =<\/span> 1<\/span>
\n
\n    dx =<\/span> (<\/span>-x + k * u)<\/span> \/ tau
\n
\n    return<\/span> dx
\n
\ndef<\/span> ode2(<\/span>t,<\/span> x)<\/span>:
\n
\n    F =<\/span> 100<\/span>
\n    k =<\/span> 10<\/span>
\n    b =<\/span> 3<\/span>
\n    m =<\/span> 1<\/span>
\n
\n    dx1 =<\/span> x[<\/span>1<\/span>]<\/span>
\n    dx2 =<\/span> (<\/span>F - k*x[<\/span>0<\/span>]<\/span> - b*x[<\/span>1<\/span>]<\/span>)<\/span>\/m
\n
\n    dx =<\/span> [<\/span>dx1,<\/span> dx2]<\/span>
\n
\n    return<\/span> dx
\n
\ndef<\/span> ode3(<\/span>t,<\/span> x)<\/span>:
\n
\n    F =<\/span> 1<\/span>
\n    k =<\/span> 10<\/span>
\n    b =<\/span> 0<\/span>
\n    m =<\/span> 1<\/span>
\n
\n    omg =<\/span> np.sqrt<\/span>(<\/span>k\/m)<\/span>
\n
\n    dx1 =<\/span> x[<\/span>1<\/span>]<\/span>
\n    dx2 =<\/span> (<\/span>F*np.sin<\/span>(<\/span>omg*t)<\/span> - k*x[<\/span>0<\/span>]<\/span> - b*x[<\/span>1<\/span>]<\/span>)<\/span>\/m
\n
\n    dx =<\/span> [<\/span>dx1,<\/span> dx2]<\/span>
\n
\n    return<\/span> dx
\n
\ndef<\/span> ode4(<\/span>t,<\/span> x)<\/span>:
\n
\n    k =<\/span> 1000<\/span>
\n    b =<\/span> 1001<\/span>
\n    m =<\/span> 1<\/span>
\n
\n    dx1 =<\/span> x[<\/span>1<\/span>]<\/span>
\n    dx2 =<\/span> (<\/span>- k*x[<\/span>0<\/span>]<\/span> - b*x[<\/span>1<\/span>]<\/span>)<\/span>\/m
\n
\n    dx =<\/span> [<\/span>dx1,<\/span> dx2]<\/span>
\n
\n    return<\/span> dx
\n
\ndef<\/span> ode_sim(<\/span>)<\/span>:
\n
\n    xt2 =<\/span> integrate.solve_ivp<\/span>(<\/span>fun=<\/span>ode2,<\/span> t_span=<\/span>[<\/span>0.0<\/span>,<\/span> 10.0<\/span>]<\/span>,<\/span>
\n                             y0=<\/span>[<\/span>0.0<\/span>,<\/span> 0.0<\/span>]<\/span>,<\/span> max_step=<\/span>0.01<\/span>,<\/span>
\n                             method=<\/span>"RK45"<\/span>)<\/span>
\n
\n    xt1 =<\/span> integrate.solve_ivp<\/span>(<\/span>fun=<\/span>ode1,<\/span> t_span=<\/span>[<\/span>0.0<\/span>,<\/span> 1.0<\/span>]<\/span>,<\/span>
\n                             y0=<\/span>[<\/span>0.0<\/span>,<\/span> 0.0<\/span>]<\/span>,<\/span> max_step=<\/span>0.01<\/span>,<\/span>
\n                             method=<\/span>"RK45"<\/span>)<\/span>
\n
\n    xt3 =<\/span> integrate.solve_ivp<\/span>(<\/span>fun=<\/span>ode3,<\/span> t_span=<\/span>[<\/span>0.0<\/span>,<\/span> 10.0<\/span>]<\/span>,<\/span>
\n                             y0=<\/span>[<\/span>0.0<\/span>,<\/span> 0.0<\/span>]<\/span>,<\/span> max_step=<\/span>0.01<\/span>,<\/span>
\n                             method=<\/span>"RK45"<\/span>)<\/span>
\n
\n    xt4 =<\/span> integrate.solve_ivp<\/span>(<\/span>fun=<\/span>ode4,<\/span> t_span=<\/span>[<\/span>0.0<\/span>,<\/span> 10.0<\/span>]<\/span>,<\/span>
\n                             y0=<\/span>[<\/span>1.0<\/span>,<\/span> 0.0<\/span>]<\/span>,<\/span> max_step=<\/span>0.1<\/span>,<\/span>
\n                             method=<\/span>"BDF"<\/span>)<\/span>
\n
\n    plt.figure<\/span>(<\/span>1<\/span>)<\/span>
\n    plt.plot<\/span>(<\/span>xt2.t<\/span>,<\/span> xt2.y<\/span>[<\/span>0<\/span>]<\/span>)<\/span>
\n    plt.title<\/span>(<\/span>"Spring-mass-dumper"<\/span>)<\/span>
\n    plt.xlabel<\/span>(<\/span>"time [t]"<\/span>)<\/span>
\n    plt.ylabel<\/span>(<\/span>"x [m]"<\/span>)<\/span>
\n    plt.grid<\/span>(<\/span>)<\/span>
\n
\n    plt.figure<\/span>(<\/span>2<\/span>)<\/span>
\n    plt.plot<\/span>(<\/span>xt1.t<\/span>,<\/span> xt1.y<\/span>[<\/span>0<\/span>]<\/span>)<\/span>
\n    plt.title<\/span>(<\/span>"ode 1"<\/span>)<\/span>
\n    plt.xlabel<\/span>(<\/span>"time [t]"<\/span>)<\/span>
\n    plt.ylabel<\/span>(<\/span>"x [m]"<\/span>)<\/span>
\n    plt.grid<\/span>(<\/span>)<\/span>
\n
\n    plt.figure<\/span>(<\/span>3<\/span>)<\/span>
\n    plt.plot<\/span>(<\/span>xt3.t<\/span>,<\/span> xt3.y<\/span>[<\/span>0<\/span>]<\/span>)<\/span>
\n    plt.title<\/span>(<\/span>"Spring-mass-resonance"<\/span>)<\/span>
\n    plt.xlabel<\/span>(<\/span>"time [t]"<\/span>)<\/span>
\n    plt.ylabel<\/span>(<\/span>"x [m]"<\/span>)<\/span>
\n    plt.grid<\/span>(<\/span>)<\/span>
\n
\n    plt.figure<\/span>(<\/span>4<\/span>)<\/span>
\n    plt.plot<\/span>(<\/span>xt4.t<\/span>,<\/span> xt4.y<\/span>[<\/span>0<\/span>]<\/span>)<\/span>
\n    plt.title<\/span>(<\/span>"Spring-mass-dumper-stiff"<\/span>)<\/span>
\n    plt.xlabel<\/span>(<\/span>"time [t]"<\/span>)<\/span>
\n    plt.ylabel<\/span>(<\/span>"x [m]"<\/span>)<\/span>
\n    plt.grid<\/span>(<\/span>)<\/span>
\n
\n    plt.show<\/span>(<\/span>)<\/span>
\n    print<\/span>(<\/span>len<\/span>(<\/span>xt4.t<\/span>)<\/span>)<\/span>
\n
\n
\nif<\/span> __name__ ==<\/span> '__main__'<\/span>:
\n    ode_sim(<\/span>)<\/span><\/div><\/td><\/tr><\/tbody><\/table><\/div>\n","protected":false},"excerpt":{"rendered":"
Dnes se pod\u00edv\u00e1me na problematiku simulace oby\u010dejn\u00e9 diferenci\u00e1ln\u00ed rovnice ODE v pythonu. Budeme vyu\u017e\u00edvat modul scipy, kter\u00fd obsahuje mnoho algortim\u016f pro numerickou matematiku. Pro \u0159e\u0161en\u00ed dynamick\u00e9ho syst\u00e9m\u016f (ten je pops\u00e1n diferenci\u00e1ln\u00edmi rovnicemi) pou\u017eijeme numerickou metodu. Rozd\u00edl mezi numerick\u00fdm a analytick\u00e9m \u0159e\u0161en\u00edm je ten, \u017ee analytick\u00e9 \u0159e\u0161en\u00ed je zcela p\u0159esn\u00e9, ale lze takto vy\u0159e\u0161it jen omezenou…<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[2],"tags":[6,4,7],"_links":{"self":[{"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/posts\/185"}],"collection":[{"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/comments?post=185"}],"version-history":[{"count":31,"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/posts\/185\/revisions"}],"predecessor-version":[{"id":254,"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/posts\/185\/revisions\/254"}],"wp:attachment":[{"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/media?parent=185"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/categories?post=185"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/mechatronicz.eu\/domains\/mechatronicz.eu\/wp-json\/wp\/v2\/tags?post=185"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}