Revision history [back]

This is a report on the updated question after the replies by rburing and tmonteil.

The answer and recommendations of rburing concerning the 1-variate example worked beautifully also for the newly proposed 2D example, as follows:

x, y = var('x,y')
def f(x,y):    
    if -2.<=x<=0.:
        return 0.
    elif -2.<=y<=0.:
        return 0.
    elif 1-x-y<=0.:
        return 0.
    elif 1.<=x<=2.:
        return 0.
    elif 1.<=y<=2.:
        return 0.
    else:
        return 0.1*exp(-0.01*((x-1/3)^2+(y-1/3)^2)/(((1-x-y)^2)*(x^2)*(y^2)))
def g(x,y):
    return x
def h(x,y):
    return y
P=parametric_plot3d([g, h, f], (x, -2, 2), (y, -2, 2), plot_points=[128,128], frame=False, color="silver", opacity=1.)
P.show(viewer='threejs')

This is solving completely the part of the problem of plotting the graph. However, the use of conditional statements yielded the graph on the whole rectangle [-2,2]^2. Although we have used x and y as barycentric coordinates here, this is still not convenient to be used as a triangular patch.

The following is a suggestion how to resolve this difficulty WITHOUT ANY USE OF CONDITIONAL STATEMENTS, by using reparametrization instead.

After a Jacobi change of variables for the triangle,

x=s(1-t), y=st

the triangle is reparametrized, with the square (s,t) \in [0,1]^2 being the new parametric domain. In the new variables, the expression in the 'else' statement above becomes

0.1exp(-0.01((1-t-1/(3.s))^2+(t-1/(3.s))^2)/(((1-s)^2)(s^2)((1-t)^2)*(t^2)))

and there follows the new code (this time including also s and t coordinate lines):

s, t = var('s,t')
def F(s,t):
    return 0.05*exp(-0.005*((1-t-1/(3.*s))^2+(t-1/(3.*s))^2)/(((1-s)^2)*(s^2)*((1-t)^2)*(t^2)))
def G(s,t):
    return s*(1-t)
def H(s,t):
    return s*t
Q0=parametric_plot3d([G, H, F], (s, 0, 1), (t, 0, 1), plot_points=[128,128], frame=False, color="gray", opacity=1.)
eps=0.001
M=64
N=64
L=128
m=var('m', domain='integer')
n=var('n', domain='integer')
curve3dplot_uline=[0 for i in range(M)] 
curve3dplot_vline=[0 for j in range(N)]
curve3dplot_uline=[parametric_plot3d([G(s,eps+m*(1.-2.*eps)/M), H(s,eps+m*(1.-2.*eps)/M), F(s,eps+m*(1.-2.*eps)/M)],  (s, 0, 1), plot_points=128, color='silver', thickness=.5, opacity=1.) for m in range(M)] 
curve3dplot_vline=[parametric_plot3d([G(eps+n*(1.-2.*eps)/N,t), H(eps+n*(1.-2.*eps)/N,t), F(eps+n*(1.-2.*eps)/N,t)],  (t, 0, 1), plot_points=128, color='silver', thickness=.5, opacity=1.) for n in range(N)]
Q1=sum(curve3dplot_uline[i] for i in range(M)) 
Q2=sum(curve3dplot_vline[j] for j in range(N))
Q=Q0+Q1+Q2
Q.show(viewer='threejs')

In this solution all conditional statements have been eliminated, the resulting graph is the graph of a parametric triangular patch. Many years ago, underflow error would have been reported at the boundary of the triangle, but numpy handles the exception well and interprets the value as 0., as needed. However, the computation of the coordinate s- and t-lines requires the introduction of a small positive parameter eps; otherwise, if eps=0, the result will be error of symbolic division by zero.

This construction by parametrization avoiding the use of conditional constraints can be straightforwardly generalized for a tetrahedron in 3D (and generally for a simplex in n-dimensions, n=2,3,4,...) . The needed generalization of the Jacobi reparametrization to n-dimensional simplex is well known (see, e.g. the 3rd volume of Fihtenholz).