def NonNegativeSolutionInRR(li):
for x in li:
if x.rhs() in RR and x.rhs()>=0:
return N(x.rhs())
|
|
BEISPIEL 4
%cython
from math import exp as fast_exp
from sage.all import srange
def discrete_sum(f,double a,double b,double step,double OberUnter): #cython anti-aging
return sum([f(c+step*OberUnter)*step for c in srange(a,b+step*OberUnter,step)])
default="math.exp(-x**2)" #monotone decreasing and symmetric
default_digs=4
def main(func=default,digs=default_digs):
def f(x):
return eval(func)
print "obersumme: ",N(discrete_sum(f,0,1,10^-digs,1)*2,digits=digs+1)
print "untersumme:",N(discrete_sum(f,0,1,10^-digs,0)*2,digits=digs+1)
main()
#@interact
def gui(func=default,digs=default_digs):
main(func,digs)
obersumme: 1.4937 untersumme: 1.4937 obersumme: 1.4937 untersumme: 1.4937 |
BEISPIEL 5
%cython
def trapez_sum(f,double n,double a,double b): #cython anti-aging
h=(b-a)/n
return h*((f(a)+f(b))/2+sum([f(a+c*h) for c in range(1,n)]))
default="sqrt(sin(1/x)+x)"
default_a=1
default_b=2
default_digs=8
def main(func=default,a=default_a,b=default_b,digs=default_digs):
accuracy=10^-digs
h=var('h'); x=var('x')
f(x)=eval(func)
MAX=f.derivative(2).find_maximum_on_interval(a,b)[1] #will find analytical solutions in future too
h=N(NonNegativeSolutionInRR(solve([accuracy==(b-a)/12*h^2*MAX],h)))
print "h needs to be smaller or equal to:",h
n=ceil(1/(h/(b-a)))
print "trapezsumme: ",N(trapez_sum(f,n,a,b),digits=digs+2)
main()
#@interact
def gui(func=default,a=default_a,b=default_b,digs=default_digs):
main(func,a,b,digs)
h needs to be smaller or equal to: 0.000318369865152595 trapezsumme: 1.458949377 h needs to be smaller or equal to: 0.000318369865152595 trapezsumme: 1.458949377 |
BEISPIEL 6
%cython
def fass_sum(f,double n,double a,double b): #cython anti-aging
h=(b-a)/n
return h/3*((f(a)+f(b))/2+sum([f(a+c*h) for c in range(n)])+2*sum([f((a+c*h+a+(c-1)*h)/2) for c in range(n+1)]))
default="sqrt(sin(1/x)+x)"
default_a=1
default_b=2
default_digs=8
def main(func=default,a=default_a,b=default_b,digs=default_digs):
accuracy=10^-digs
h=var('h'); x=var('x')
f(x)=eval(func)
MAX=f.derivative(4).find_maximum_on_interval(a,b)[1] #will find analytical solutions in future too
h=N(NonNegativeSolutionInRR(solve([accuracy==(b-a)/2880*h^4*MAX],h)))
print "h needs to be smaller or equal to:",h
n=ceil(1/(h/(b-a)))
print "fasssumme: ",N(fass_sum(f,n,a,b),digits=digs+2)
main()
#@interact
def gui(func=default,a=default_a,b=default_b,digs=default_digs):
main(func,a,b,digs)
h needs to be smaller or equal to: 0.0639970891248658 fasssumme: 1.543543001 h needs to be smaller or equal to: 0.0639970891248658 fasssumme: 1.543543001 |
BEISPIEL 7
%cython
from sage.all import srange
def ellipse_approx(double step,double area,f): #cython anti-aging
lastY=None; len=0; n=2
for x in srange(-area+step/2,area,step):
Y=f(x)
if Y!=None:
if lastY!=None:
len+=sqrt(step**2+(Y-lastY)**2)
n+=2
else:
len+=Y
lastY=Y;
len+=lastY
return 2*len,n
default="sqrt((-2*x**2+5)/3)"
default_steps=0.5
default_area=2
default_draw=true
def main(glt=default,step=default_steps,area=default_area,draw=default_draw):
x=var('x'); y=var('y'); gl=eval(glt)
def liY(x):
try:
return eval(glt)
except:
return None
re=ellipse_approx(step,area,liY)
print "umfang",N(re[0])
print "line segments:",re[1]
if draw:
g=implicit_plot(eval("y**2==("+glt+")*("+glt+")"),(x,-area,area),(y,-area,area))
for x in srange(-area+step/2,area,step):
g+=line([(x,-area),(x,area)])
show(g,aspect_ratio=1)
main()
#@interact
def gui(glt=default,step=default_steps,area=default_area,draw=default_draw):
main(glt,step,area,draw)
umfang 8.66935857497510 line segments: 12 umfang 8.66935857497510 line segments: 12 |
