Ticket #10983: graphique_3_fixed.py

File graphique_3_fixed.py, 4.4 KB (added by casamayou, 5 years ago)
Line 
1r"""
2sage: import scipy; from scipy import integrate
3sage: f = lambda y, t: - cos(y * t)
4sage: t = srange(0, 5, 0.1); p = Graphics()
5sage: for k in srange(0, 10, 0.15):
6...      y = integrate.odeint(f, k, t)
7...      p += line(zip(t, flatten(y)))
8sage: t = srange(0, -5, -0.1); q = Graphics()
9sage: for k in srange(0, 10, 0.15):
10...      y = integrate.odeint(f, k, t)
11...      q += line(zip(t, flatten(y)))
12sage: y = var('y')
13sage: v = plot_vector_field((1, -cos(x * y)), (x,-5,5), (y,-2,11))
14sage: g = p + q + v; # g.show()
15"""
16
17r"""
18sage: import scipy; from scipy import integrate
19sage: a, b, c, d = 1., 0.1, 1.5, 0.75
20sage: def dX_dt(X, t=0):
21...    return [ a*X[0] -   b*X[0]*X[1] ,
22...            -c*X[1] + d*b*X[0]*X[1] ]
23sage: t = srange(0, 15, .01)                     # echelle de temps
24sage: X0 = [10, 5]  # conditions initiales : 10 lapins et 5 renards
25sage: X = integrate.odeint(dX_dt, X0, t)     # resolution numerique
26sage: lapins, renards =  X.T         # raccourcis de  X.transpose()
27sage: p = line(zip(t, lapins), color='red') # trace du nb de lapins
28sage: p += text("Lapins",(12,37), fontsize=10, color='red')
29sage: p += line(zip(t, renards), color='blue')# idem pr les renards
30sage: p += text("Renards",(12,7), fontsize=10, color='blue')
31sage: p.axes_labels(["temps", "population"]); # p.show(gridlines=True)
32sage: ### Deuxieme graphique :
33sage: n = 11;  L = srange(6, 18, 12 / n); R=srange(3, 9, 6 / n)
34sage: CI = zip(L, R)                # liste des conditions initiales
35sage: def g(x,y):
36...       v = vector(dX_dt([x, y]))  # pour un trace plus lisible,
37...       return v/v.norm()          # on norme le champ de vecteurs
38sage: x, y = var('x, y')
39sage: q = plot_vector_field(g(x, y), (x, 0, 60), (y, 0, 36))
40sage: for j in range(n):
41...       X = integrate.odeint(dX_dt, CI[j], t)        # resolution
42...       q += line(X, color=hue(.8-float(j)/(1.8*n))) #  graphique
43sage: q.axes_labels(["lapins","renards"]); # q.show()
44"""
45
46r"""
47sage: x, y, t = var('x, y, t')
48sage: alpha(t) = 1; beta(t) = t/2; gamma(t) = t + t**3/8
49sage: env = solve([alpha(t)*x+beta(t)*y==gamma(t),\
50...         diff(alpha(t),t)*x+diff(beta(t),t)*y==diff(gamma(t),t)],\
51...         [x,y])
52sage: f = lambda x:x^2 / 4
53sage: p = plot(f, -8, 8, rgbcolor=(0.2,0.2,0.4)) # trace la parabole
54sage: for u in srange(0, 8, 0.1): # trace des normales a la parabole
55...       p += line([[u, f(u)], [-8*u, f(u) + 18]], thickness=.3)
56...       p += line([[-u, f(u)], [8*u, f(u) + 18]], thickness=.3)
57sage: p += parametric_plot((env[0][0].rhs(),env[0][1].rhs()),\
58...       (t, -8, 8),color='red')              # trace la developpee
59sage: # p.show(xmin=-8, xmax=8, ymin=-1, ymax=12, aspect_ratio=1)
60sage: t = var('t'); p = 2
61sage: x(t) = t; y(t) = t^2 / (2 * p)
62sage: f(t) = [x(t), y(t)]
63sage: df(t) = [x(t).diff(t), y(t).diff(t)]
64sage: d2f(t) = [x(t).diff(t, 2), y(t).diff(t, 2)]
65sage: T(t) = [df(t)[0] / df(t).norm(), df[1](t) / df(t).norm()]
66sage: N(t) = [-df(t)[1] / df(t).norm(), df[0](t) / df(t).norm()]
67sage: R(t) =  (df(t).norm())^3 / \
68...        (df(t)[0] * d2f(t)[1] -df(t)[1] * d2f(t)[0])
69sage: Omega(t) = [f(t)[0] + R(t)*N(t)[0], f(t)[1] + R(t)*N(t)[1]]
70sage: g = parametric_plot(f(t), (t, -8, 8), color='green', thickness=2)
71sage: for u in srange(.4, 4, .2):
72...    g += line([f(t = u), Omega(t = u)], color='red', alpha = .5)
73...    g += circle(Omega(t = u), R(t = u), color='blue')
74sage: # g.show(aspect_ratio=1, xmin=-12, xmax=7, ymin=-3, ymax=12)
75"""
76
77r"""
78sage: u, v = var('u, v')
79sage: h = lambda u,v: u^2 + 2*v^2
80sage: f = plot3d(h, (u,-1,1), (v,-1,1), aspect_ratio=[1,1,1])
81sage: f(x, y) = x^2 * y / (x^4 + y^2)
82sage: t, theta = var('t theta')
83sage: limit(f(t * cos(theta), t * sin(theta)) / t, t=0)
84cos(theta)^2/sin(theta)
85sage: solve(f(x,y) == 1/2, y)
86[y == x^2]
87sage: a = var('a'); h = f(x, a*x^2).simplify_rational(); h
88a/(a^2 + 1)
89sage: g = plot(h, a, -4, 4)
90sage: p = plot3d(f(x, y), (x,-2,2), (y,-2,2), plot_points=[150,150])
91sage: for i in range(1,4):
92...    p += plot3d(-0.5 + i / 4, (x, -2, 2), (y, -2, 2),\
93...                 color=hue(i / 10), opacity=.1)
94"""
95
96r"""
97sage: x, y, z = var('x, y, z'); a = 1
98sage: h = lambda x, y, z:(a^2 + x^2 + y^2)^2 - 4*a^2*x^2-z^4
99sage: f = implicit_plot3d(h, (x, -3, 3), (y, -3, 3), (z, -2, 2),\
100...                    plot_points=100, adaptative=True)
101sage: g1 = line3d([(-10*cos(t)-2*cos(5*t)+15*sin(2*t),\
102...     -15*cos(2*t)+10*sin(t)-2*sin(5*t),\
103...      10*cos(3*t)) for t in srange(0,6.4,.1)],radius=.5)
104
105"""