I create a list Q (below) of pairs of numbers [t1,z1], I want the pairs (maybe just a print of them?) at the local extrema (of z1) and where z1 crosses zero. Is there a simple way to do this?
var('t y z')
P=desolve_system_rk4([z,0.005*exp(-y/8000)*z^2-10],[y,z],ics=[0,32000,0],ivar=t,end_points=400)
Q=[ [t1,z1] for t1,y1,z1 in P]
Thanks, Wayne
The code:
var('t y z')
P = desolve_system_rk4( [ z, 0.005 * exp(-y/8000) * z^2 - 10]
, [y,z]
, ics = [ 0, 32000, 0 ]
, ivar=t
, end_points=400 )
Q = [ (t1,z1) for t1,y1,z1 in P ]
for k in xrange( 1, len(Q)-1 ):
k1, k2, k3 = k-1, k, k+1
(t1,z1), (t2,z2), (t3,z3) = Q[k1], Q[k2], Q[k3]
if z1 <= z2 and z2 >= z3:
print "( %s, %s ) is a local MAX" % ( t2, z2 )
if z1 >= z2 and z2 <= z3:
print "( %s, %s ) is a local MIN" % ( t2, z2 )
if z1 <= 0 and z2 >=0:
print "ZERO BETWEEN ( %s, %s ) and ( %s, %s )" % ( t1, z1, t2, z2 )
if z1 >= 0 and z2 <=0:
print "ZERO BETWEEN ( %s, %s ) and ( %s, %s )" % ( t1, z1, t2, z2 ) # EDITED HERE, sorry
finds...
ZERO BETWEEN ( 0, 0 ) and ( 0.1, -0.9999969473920776 )
( 37.8, -231.6456047650853 ) is a local MIN
but misses of course a possible zero at the end of the interval. Bur looking at
sage: Q[-2:]
[(399.9, -36.4642141232692), (400.0, -36.4558355793491)]
we see there is no zero at the other end of the interval.
Another option is to use splines to approximate the curve. Here is some code inspired by https://stackoverflow.com/questions/3...
var('t y z')
P=desolve_system_rk4([z,-z^2*(-2+t/(1+t))-y],[y,z],ics=[0.1,10,2],ivar=t,end_points=100)
Q=[ [float(t1),float(z1)] for t1,y1,z1 in P]
x_points=([p[0] for p in Q])
y_points=([p[1] for p in Q])
tmp = interpolate.splrep(x_points, y_points)
f=lambda x: interpolate.splev(x, tmp)
You could then find a zero using:
find_root(f,12,20)
Also, you can compute the derivative using:
g=lambda x: interpolate.splev(x, tmp,der=1)
Then, find potential extrema using find_root.
OK, I think I have a path to a solution. I have changed the differential equation so that it has zeros and this code prints the values on either side of where the function crosses zero. I can do something similar for the extrema. So I guess I am now asking if there is a function call that does this, or if there is a smarter way. Thanks, Wayne
var('t y z')
P=desolve_system_rk4([z,-z^2*(-2+t/(1+t))-y],[y,z],ics=[0.1,10,2],ivar=t,end_points=100)
Q=[ [t1,y1] for t1,y1,z1 in P]
line(Q).show()
print len(Q)
[a1,b1]=Q[0]
for i in range(1,len(Q)):
[a2,b2]=Q[i]
if(((b1<0.)and(b2>0.)) or ((b1>0.)and(b2<0.))):
print a1,b1,a2,b2
a1=a2
b1=b2
Asked: 2018-01-11 13:27:25 +0200
Seen: 480 times
Last updated: Jan 13 '18
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.