# Ticket #10983: graphique_3.py

File graphique_3.py, 4.4 KB (added by casamayou, 9 years ago)

updated file

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