解微分方程组,变某一个参数

par_range.py

 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
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import matplotlib as mpl
mpl.rcParams["font.family"] = "Helvetica Neue"

from par_range_sm import init, par, rep, par_rand

t = np.arange(0, 30, 0.01)
init = init()

for i in range(4):
    a1 = par_rand(100)
    param = par(a1=a1)
    track = odeint(rep, init, t, args=(param,))
    M1, M2, M3, P1, P2, P3 = track.T
    filename = "a1_%05.1f" % a1
    pdf = filename + ".pdf"
    png = filename + ".png"

    plt.figure(figsize=(8, 6), dpi=80)
    plt.title("$a_1=%5.1f$" % a1, x=0.4, y=1.02, ha="left")
    plt.axis([0, 30, 0, 60])
    plt.xlabel("Time")
    plt.ylabel("Concentration")
    plt.xticks([0, 10, 20, 30])
    plt.yticks([0, 20, 40, 60])
    plt.plot(t, M1, "black", lw=2, ls="-", label="$M_1$")
    plt.plot(t, M2, "red", lw=2, ls="--", label="$M_2$")
    plt.plot(t, M3, "blue", lw=2, ls="-.", label="$M_3$")
    plt.legend(loc=9, ncol=3)
    plt.savefig(png)
    plt.close()
    print "%s saved!" % png

print "Done!"

其中 init, param, reppar_rand 函数来自另一个文件 par_range_sm.py

 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
import numpy as np


def init(m1=40,
         m2=2,
         m3=2,
         p1=0,
         p2=2,
         p3=2):
    return (m1, m2, m3, p1, p2, p3)


def par(a1=100,
        a2=100,
        a3=100,
        b1=5,
        b2=5,
        b3=5,
        r1=5,
        r2=5,
        r3=5,
        m=2):
    return (a1, a2, a3, b1, b2, b3, r1, r2, r3, m)


def rep(init, t, param):
    (m1, m2, m3, p1, p2, p3) = init
    (a1, a2, a3, b1, b2, b3, r1, r2, r3, m) = param

    M1 = -m1 + a1/(1+p3**m)
    M2 = -m2 + a2/(1+p1**m)
    M3 = -m3 + a3/(1+p2**m)
    P1 = b1*m1 - r1*p1
    P2 = b2*m2 - r2*p2
    P3 = b3*m3 - r3*p3
    return np.array([M1, M2, M3, P1, P2, P3])


def par_rand(a, sigma=0.15):
    return np.random.uniform(
        a * (1-sigma),
        a * (1+sigma)
        )

示例:

Image Title

2014-05-09 16:0046PythonScipyNumpyMatplotlib
comments powered by Disqus