florin's profile - activity

Thanks for the fixup of parametricplot ! indeed the three branches appear only after enlarging the t range to say -8,8,

Thanks for the fixup of parametricplot ! indeed the three branches appear only after enlarging the t range to say -8,8,

Thanks for the fixup of parametricplot ! indeed the three branches appear only after enlarging the t range to say -8,8,

Hi The following plot

var('t')
parametric_plot( (8*(6+t)/(t^2-16), 2*(8+ 3 *t)/(t^2-16)), (t, -6, 6),xmin=-3,xmax=3,detect_poles=True)

has three branches. 1) For some reason, the plot is not correct. Also

var('t')
parametric_plot( (8*(6+t)/(t^2-16), 2*(8+ 3 *t)/(t^2-16)), (t, -6, 6),xmin=-3,xmax=3,ymin=-3,ymax=3)

fails. 2)How to plot the branches in different colors? 3) Can I add a point moved by a cursor when t increases, like with the Mathematica command Manipulate? Or make an animation to show how a point traverses the three branches?

A parametric plot with branches for which restricting the range with xmin, xmax doesn't work Hi The following plot v

A parametric plot with branches for which restricting the range with xmin, xmax doesn't work Hi The following plot v

How to restrict the range of a parametric plot and draw the branches in different colors Hi The following plot var('

How to restrict the range of a parametric plot and draw the branches in different colors Hi The following plot var('

How to restrict the range of a parametric plot and draw the branches in different colors Hi The following plot var('

draw a parametric plot with several branches in different colors Hi The following plot var('t') parametric_plot( (8*(6

transportation problem asks finding the least cost solution for moving merchandize between "supplies" and demands". It

Merci, Emmanuel, now it works :)

def optL(f,h):
   var('a, b, x, y, lam')
   assume(a>=0,b>=0,x>= 0)
   L = f + h * lam
   DL = L.gradient([x, y])
   sol = solve([DL[k] == 0 for k in range(len(DL))]+[ h == 0],x, y, lam,
   solution_dict=True)
   ti=[["$\\lambda$","$x$", "$y$", "$f(x, y)$"]]
   re=[(sol[j][lam],sol[j][x], sol[j][y],f(x=sol[j][x], y=sol[j][y])) for j in range(len(sol))]
   ans=table(ti + re, header_row = True)
   show(ans)
   return re
ex5=optL(x*y,x^2/a^2+ y^2/b^2-1)
ex5

Hi I have a sequence of commands which does Lagrangian optimization

var('a, b, x, y, lam')
assume(a>=0,b>=0,x>= 0)
f = x*y #minimiser et maximiser sous contrainte d'egalite h=0
h = x^2/a^2 + y^2/b^2-1
L = f + h * lam
DL = L.gradient([x, y])
sol = solve([DL[k] == 0 for k in range(len(DL))]+[ h == 0],x, y, lam,
solution_dict=True)
ti=[["$\\lambda$","$x$", "$y$", "$f(x, y)$"]]
re=[(sol[j]
[lam],sol[j][x], sol[j][y],f(x=sol[j][x], y=sol[j][y])) for j in
range(4)]
table(ti + re, header_row = True)

and I would like to turn it into a function. I add the declaration, but I get no output. I know neither why, nor how to debug this. Thanks for any help :)

def optL(f,h):
       var('a, b, x, y, lam')
       assume(a>=0,b>=0,x>= 0)
       L = f + h * lam
       DL = L.gradient([x, y])
       sol = solve([DL[k] == 0 for k in range(len(DL))]+[ h == 0],x, y, lam,
       solution_dict=True)
       ti=[["$\\lambda$","$x$", "$y$", "$f(x, y)$"]]
       re=[(sol[j][lam],sol[j][x], sol[j][y],f(x=sol[j][x], y=sol[j][y])) for j in range(4)]
       table(ti + re, header_row = True)
optL(x*y,x^2/a^2 + y^2/b^2-1)

The answer is InteractiveLPProblem??

I'm using sage 8.8 on Windows 10. Thanks !

I would like something using the special structure of the problem, for teaching purposes, like the wonderful InteractiveLPProblem for linear programming. Thanks

Hi I tried both html and .pdf, but the 3dplots seem to have disappeared. Is there a trick to get a copy of your session, includimg the 3dplots?

Thanks! The error on 8.8 came maybe from Spectral

I reinstalled the latest Windows version 8.8. repeating now a second time the show(g1) yields error!

RuntimeError                              Traceback (most recent call last)
<ipython-input-20-2dbcc4c9974d> in <module>()
----> 1 show(g1)

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/repl/rich_output/pretty_print.pyc in show(*args, **kwds)
    256         args[0].show()
    257         return
--> 258     pretty_print(*args, **kwds)

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/repl/rich_output/pretty_print.pyc in pretty_print(*args, **kwds)
    227             pass
    228         elif len(args) == 1:
--> 229             dm.display_immediately(*args, **kwds)
    230         else:
    231             SequencePrettyPrinter(*args, **kwds).pretty_print()

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/repl/rich_output/display_manager.pyc in display_immediately(self, obj, **rich_repr_kwds)
    833             1/2
    834         """
--> 835         plain_text, rich_output = self._rich_output_formatter(obj, rich_repr_kwds)
    836         self._backend.display_immediately(plain_text, rich_output)
    837 

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/repl/rich_output/display_manager.pyc in _rich_output_formatter(self, obj, rich_repr_kwds)
    633         if rich_output is None:
    634             rich_output = self._preferred_text_formatter(
--> 635                 obj, plain_text=plain_text, **rich_repr_kwds)
    636         # promote output container types to backend-specific containers
    637         plain_text = self._promote_output(plain_text)

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/repl/rich_output/display_manager.pyc in _preferred_text_formatter(self, obj, plain_text, **kwds)
    535             return out
    536         if want == 'latex' and OutputLatex in supported:
--> 537             out = self._backend.latex_formatter(obj, **kwds)
    538             if type(out) is not OutputLatex:
    539                 raise OutputTypeException('backend returned wrong output type, require Latex')

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/repl/rich_output/backend_base.pyc in latex_formatter(self, obj, **kwds)
    482             mathjax = MathJax().eval(obj, mode='plain', combine_all=True)
    483         else:
--> 484             mathjax = MathJax().eval(obj, mode='plain', combine_all=False)
    485         from sage.repl.rich_output.output_basic import OutputLatex
    486         return OutputLatex(str(mathjax))

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/misc/latex.pyc in eval(self, x, globals, locals, mode, combine_all)
   1951         """
   1952         # Get a regular LaTeX representation of x
-> 1953         x = latex(x, combine_all=combine_all)
   1954 
   1955         # The following block, hopefully, can be removed in some future MathJax.

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/misc/latex.pyc in __call__(self, x, combine_all)
    923         """
    924         if has_latex_attr(x):
--> 925             return LatexExpr(x._latex_())
    926         try:
    927             f = latex_table[type(x)]

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/plot/graphics.py in _latex_(self, **kwds)
   3240         """
   3241         tmpfilename = tmp_filename(ext='.pgf')
-> 3242         self.save(filename=tmpfilename, **kwds)
   3243         with open(tmpfilename, "r") as tmpfile:
   3244                 latex_list = tmpfile.readlines()

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/misc/decorators.pyc in wrapper(*args, **kwds)
    410             kwds[self.name + "options"] = suboptions
    411 
--> 412             return func(*args, **kwds)
    413 
    414         # Add the options specified by @options to the signature of the wrapped

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/plot/graphics.py in save(self, filename, **kwds)
   3165             rc_backup = (rcParams['ps.useafm'], rcParams['pdf.use14corefonts'],
   3166                          rcParams['text.usetex']) # save the rcParams
-> 3167             figure = self.matplotlib(**options)
   3168             # You can output in PNG, PS, EPS, PDF, PGF, or SVG format, depending
   3169             # on the file extension.

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/plot/graphics.py in matplotlib(self, filename, xmin, xmax, ymin, ymax, figsize, figure, sub, axes, axes_labels, axes_labels_size, fontsize, frame, verify, aspect_ratio, gridlines, gridlinesstyle, vgridlinesstyle, hgridlinesstyle, show_legend, legend_options, axes_pad, ticks_integer, tick_formatter, ticks, title, title_pos, base, scale, stylesheet, typeset)
   2616                     pass
   2617             g.set_options(opts)
-> 2618             g._render_on_subplot(subplot)
   2619             if hasattr(g, '_bbox_extra_artists'):
   2620                 self._bbox_extra_artists.extend(g._bbox_extra_artists)

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/plot/contour_plot.py in _render_on_subplot(self, subplot)
    159         contours = options['contours']
    160         if 'cmap' in options:
--> 161             cmap = get_cmap(options['cmap'])
    162         elif fill or contours is None:
    163             cmap = get_cmap('gray')

/opt/sagemath-8.8/local/lib/python2.7/site-packages/sage/plot/colors.py in get_cmap(cmap)
   1418     elif isinstance(cmap, six.string_types):
   1419         if not cmap in cm.datad:
-> 1420             raise RuntimeError("Color map %s not known (type import matplotlib.cm; matplotlib.cm.datad.keys() for valid names)" % cmap)
   1421         return cm.__dict__[cmap]
   1422 

RuntimeError: Color map spectral not known (type import matplotlib.cm; matplotlib.cm.datad.keys() for valid names)

Hi The code

f(x,y)=x^3-x^2 + y^2;print f
g1 = contour_plot(f(x,y),(x,-2,2),(y,-2,2),contours=20,fill=False,cmap='spectral')
show(g1)

produces a plot with a hyperbolic point. On my computer, repeating this

show(g1)

yields a different plot

The purpose was to show this together with gradient field

f(x,y)=(x+2)^2/(y^2+ (x-1/2)^2) ;print f
g1 = contour_plot(f(x,y),(x,-2,6),(y,-2,4),contours=[1/10,1/2,4/5,1,2,3],cmap='spectral',fill=False)
df=f.gradient()
show(df)
g2 = plot_vector_field(df,(x,-2,2),(y,-2,2))

This doesn't work, since the first plot is incorrect show(g1+g2)

Hi, I just spent half an hour on these two commands.

def mov(): 
        P=matrix(QQ,[[1, 1],[0, 1]]);
        show(P)

I retyped them, I moved them in different cells, nothing, I get syntax error.

"<ipython-input-76-049764d5c40b>", line 1
    def mov() :
       ^
SyntaxError: invalid syntax

It works as long as I do not try to make a function

Sorry, I'm missing some syntax. I tried

def LP():
    alias 'LP'='MixedIntegerLinearProgram'
    load("LP.sage")

gets the error

File "<ipython-input-4-fc93b025772b>", line 2
    alias 'LP'='MixedIntegerLinearProgram'
             ^
SyntaxError: invalid syntax

Hi. This morning I had two similar codes that worked:

var('s')
L=(1-exp(-s))/(s)
Ls=L.simplify_full()
show(parent(Ls))
tay=Ls.series(s,2*n+1) #or tay=taylor(Ls,s,0,2*n+1)?

and

    var('s')
    L=(1-exp(-s))/(s)
    Ls=L.simplify_full()
    show(parent(Ls))
    tay=taylor(Ls,s,0,2*n+1)

I couldn't see any difference in the parents, but with series I could continue getting Pade approx:

ta=tay.power_series(QQbar)
show(parent(ta))
n=2;Lp=ta.pade(n-1,n)

and with taylor I couldn't, hence my first question. Now, not even the first part, getting the series, does not work anymore, which adds a second question. It is true that the code that works was more complicated (including extra manipulations :the Pollaczek-Khinchine formula)

var('s ');rho=2/3;
#Compute Pollackek-Khinchine formula L_rui for the Laplace transf of ruin prob
R.<s> = PolynomialRing(QQbar) #algebraic numbers over Q, (or RR,CC)
FF = R.fraction_field()
L_f=(1-exp(-s))/(s);L_F=(1-L_f)/s
#Compute Pollackek-Khinchine formula L_rui for the Laplace transf of ruin prob
m1=limit(L_F,s=0) #the s is removable singularity
fe=L_F/m1
Fe=(1-fe)/s
L_rui=rho*Fe/(1-rho*fe) #in SR
Ls=L_rui.simplify_full()
n=2;
show(taylor(Ls,s,0, 2*n+1))
tay=Ls.series(s,2*n+1)
show(parent(tay))
ta=tay.power_series(QQbar)
show(parent(ta))
Lp=ta.pade(n-1,n)
wh,dec=Lp.partial_fraction_decomposition()
dec

Clearly, when simplifying for the purpose of asking on this forum, I added some silly mistake, but I can't see which

The working code for n=3 is

A=[1/6,2/6,3/6];ex=[2,6,4];rho=5/8;var('x s u')
F=sum(A[i]*exp(-x*ex[i]) for i in [0..2])#sum of exponentials
R.<s> = PolynomialRing(QQbar)
FF = R.fraction_field()
L_F=sum(A[i] /(s+ex[i]) for i in [0..2]);#Laplace transform of F
#Compute Pollackek-Khinchine formula L_rui for the Laplace transform of ruin probability
m1=L_F(s=0)
fe=L_F/m1
Fe=(1-fe)/s 
show(Fe) #note the s is simplified when working in QQbar
L_rui=rho*Fe/(1-rho*fe)
R.change_ring(RR)
FF = R.fraction_field()
wh,LL=L_rui.partial_fraction_decomposition()
show(LL[0])

with output

[0.5332274783298694?s+1.087668109036214?,0.07804278965136865?s+2.967947295921641?,0.01372973201876195?s+5.575963542410567?]

I need now to "remove" somehow the ? before taking

inverse_laplace(LL[0],s,u)

It would be nice if something like

inverse_laplace(LL[0].change_ring(RR),s,u)

or inverse_laplace(RR(LL[0]),s,u) would work

Hi, I'm trying to invert the "Pollaczek-Khinchine" Laplace transform when it is rational

This works for me at degree 2:

var('x,s')
Fx = (1/6*exp(-2*x)+5/6*exp(-6*x));rho=2/3
print('Hyperexponential claims:',Fx)
R.<s> = PolynomialRing(QQbar)#when all coefficients are not integer, use CC
FF = R.fraction_field()
L_F=laplace(Fx,x,s)#Laplace transform of F
#Compute Pollackek-Khinchine (PK) formula L_rui for the Laplace transform (LT) of ruin probability
m1=L_F(s=0)
fe=L_F/m1
Fe=(1-fe)/s
L_rui=rho*Fe/(1-rho*fe)
show(L_rui.simplify_full())
inverse_laplace(SR(L_rui),s,u)

but not at degree 3, since I do not know how to use partial_fraction_decomposition, and then to switch to RR numbers and then invert . If I start in R. = PolynomialRing(RR), for an already known LT, everything is fine. But, a certain simplification by s in PK formula will become impossible due to rounding errors, so I am forced to start with R. = PolynomialRing(QQbar) After obtaining the partial_fraction_decomposition, I must apply RR to all numbers , but I do not manage to do it. Without that conversion, inverse_laplace won't work

OK I got it, I just need to replace QQbar by RR. Silly me !!!!

the conversion is doable by RR(). Is there an easy way to apply RR to all parts of an expression?

The Pade of the uniform program is

var('s'); rho=2/3;L_f=(1-exp(-s))/(s) #LT of uniform density
L_F=(1-L_f)/s #LT of uniform survival function
#Compute Pollackek-Khinchine formula L_rui for the Laplace transform of ruin probability
m1=limit(L_F,s=0) #the s is removable singularity
fe=L_F/m1
Fe=(1-fe)/s 
L_rui=rho*Fe/(1-rho*fe)
show(L_rui)
n=2;
tay=taylor(L_rui,s,0,2*n+1)
ta=tay.power_series(QQ)
Lp=ta.pade(1,2)
Lp

with answer

(5/8*s + 15/2)/(s^2 + 45/4*s + 45/4)

it only remains to parfrac this. Thanks, Florin

in competitive risks models, there are more interesting things to do OK, I'm dealing with a risk expert :) which University are you based in?

I am just trying now to program in Sage a pedagogic exercise: compute the ruin probability for the Cramer-Lundberg model with a) hyperexponential claims; this should be easy, except that I am novice enough to Sage to have forgotten how these 1.33456? numbers are called, and how to convert an expression involving many of these numbers to a type which will be accepted by the command inverse_laplace

b) uniform claims (so the LT is not rational); after rationalizing the LT by Pade, this reduces to previous; it's a one liner in Mathematica, but takes hours in sage due to the crude help system.