In [1]:
%matplotlib inline
In [2]:
#import math
#import numpy
#from numpy import *
#import scipy
from scipy import integrate
from pylab import *
In [3]:
y = lambda x: x*sin(x)
print y(2)
In [7]:
integrate.quad(y,0,3)
Out[7]:
In [7]:
print integrate.quad(y,0,5)
print integrate.quad(y,0,20)
In [8]:
print(-5*cos(5)+sin(5))
In [9]:
integrate.odeint?
In [14]:
def deriv(y,t): # return derivatives of the array y
a = -4.0
b = -0.1
return array([ y[1], a*y[0]+b*y[1] ])
time = linspace(0.0,30.0,1000)
yinit = array([0.0005,0.2]) # initial values
y = integrate.odeint(deriv,yinit,time)
figure()
plot(time,y[:,0],lw=2) # y[:,0] is the first column of y
plot(time,y[:,1],lw=2)
#xlabel(‘t’)
#ylabel(‘y’)
show()
In [15]:
def SIR(y,t):
r = 0.001
rho = 0.9
return array([ -r*y[0]*y[1], r*y[0]*y[1] - rho*y[1], rho*y[1]])
time1 = linspace(0,5,100)
time2 = linspace(0,5,1000)
yinit = array([ 9000, 1000, 0]) # Initial values of S, I and R.
y2 = integrate.odeint(SIR, yinit, time2)
plot(time2, y2[:,0], lw=2)
plot(time2, y2[:,1], lw=2)
plot(time2, y2[:,2], lw=2)
xlabel("time")
ylabel("number of people")
legend(['susceptible', 'infected', 'recovered'], loc='best')
Out[15]:
In [16]:
k = 0.7
y0 = 1.0
def f(x,t,k):
return k*x
t = linspace(0,2,1000)
y1 = integrate.odeint(f, y0, t, (k,))
y2 = integrate.odeint(f, 8, t, (-0.5,))
figure()
plot(t, y1, lw=2)
plot(t, y2, lw=2)
Out[16]:
In [17]:
# time steps: an array of values starting from 0 going up to (but
# excluding) 10, in steps of 0.01
t = arange(0, 6., 0.01)
# parameters
r = 2.
K = 10.
# initial condition
x0 = 0.1
# let's define the right-hand side of the differential equation
# It must be a function of the dependent variable (x) and of the
# time (t), even if time does not appear explicitly
# this is how you define a function:
def f(x, t, r, K):
# in python, there are no curling braces '{}' to start or
# end a function, nor any special keyword: the block is defined
# by leading spaces (usually 4)
# arithmetic is done the same as in other languages: + - * /
return r*x*(1-x/K)
# call the function that performs the integration
# the order of the arguments is as below: the derivative function,
# the initial condition, the points where we want the solution, and
# a list of parameters
x = integrate.odeint(f, x0, t, (r, K))
# plot the solution
plot(t, x, lw=2)
xlabel('t') # define label of x-axis
ylabel('x') # and of y-axis
# plot analytical solution
# notice that `t` is an array: when you do any arithmetical operation
# with an array, it is the same as doing it for each element
plot(t, K * x0 * exp(r*t)/(K+x0*(exp(r*t)-1.)),lw=6,alpha=0.4)
legend(['approximation', 'analytical solution'], loc='best') # draw legend
Out[17]:
In [19]:
def f(x):
if x >= 0:
return 2.0*x**(2.0/3) + x**(5.0/3)
else:
return 2.0*(-x)**(2.0/3) - (-x)**(5.0/3)
#print f(-0.4)
t = linspace(-2,2,1000)
ft = map(f,t)
ylim(-0.5,3)
xlim(-2,1)
xlabel('x')
ylabel('y')
title('plot')
plot(t,ft,lw=3)
Out[19]:
In [113]:
def
In [13]:
import numpy as np
x = np.random.random(50)
y = np.random.random(50)
c = np.random.random(50) # color of points
s = 500 * np.random.random(50) # size of points
fig, ax = plt.subplots()
im = ax.scatter(x, y, c=c, s=s, cmap=plt.cm.jet)
# Add a colorbar
fig.colorbar(im, ax=ax)
# set the color limits - not necessary here, but good to know how.
im.set_clim(0.0, 1.0)
In [63]:
def g(b,e,f,V0,t,Q):
return b*e - f*Q/(V0 + (e-f)*t)
def A(t,a,b,e,f,V0):
return b*(V0 + (e-f)*t) + (a - b*V0)*V0**(f/(e-f))*(V0 + (e-f)*t)**(-f/(e-f))
Q0 = 100
t = linspace(3,4.5,1000)
a = A(t,10,2,3,100,10)
Q = integrate.odeint(g, Q0, t, (10,2,3,100))
plot(t,Q)
print a[2.0]
In [196]:
def decay_chain(y, t):
k1 = 0.15
k2 = 0.2
return array([-k1*y[0], k1*y[0]-k2*y[1], k2*y[1]])
t = linspace(0,20,1000)
init_vals = array([200,100,0])
y = integrate.odeint(decay_chain, init_vals, t)
plot(t, y[:,0], lw=2)
plot(t, y[:,1], lw=2)
plot(t, y[:,2], lw=2)
Out[196]:
In [69]:
def decay_chain(y, t, k1, k2):
return array([-k1*y[0], k1*y[0]-k2*y[1], k2*y[1]])
t = linspace(0,50,1000)
init_vals = array([200,100,0])
k1 = [0.05, 0.1, 0.15]
k2 = [0.05, 0.1, 0.15]
y = []
for i in range(len(k1)):
for j in range(len(k2)):
y.append(integrate.odeint(decay_chain, init_vals, t, (k1[i],k2[j])))
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(10,10))
for i in range(len(k1)):
for j in range(len(k2)):
axes[i][j].plot(t, y[3*i+j][:,0], lw=2)
axes[i][j].plot(t, y[3*i+j][:,1], lw=2)
axes[i][j].plot(t, y[3*i+j][:,2], lw=2)
axes[i][j].set_title(r"$k_1 = $" + " " + str(k1[i]) + ", "
+ r"$k_2 = $" + " " + str(k2[j]))
legend(['lead','uranium','plutonium'],loc="best",prop={"size":10})
Out[69]:
In [45]:
y = []
k1 = [0.05, 0.1, 0.15]
print len(k1)
print range(len(k1))
for i in range(len(k1)):
y.append(k1[i])
print y
In [72]:
def f(x,y):
return x*log(x**2) + log(y)
t = linspace(0.1,10,100)
plot(t,sin(t)*log(t),lw=3,color='red',label=)
Out[72]:
In []: