Advanced Search

Question: plot3d of a difficult integral

Question:plot3d of a difficult integral

GGEZ 10 Maple
plot3d
July 23 2013
1

I'm using maple 16 to do a surface plot w/ plot3d of a fairly complex integral. The integral involves a number of nested functions. The plot ends up taking a very long time. The surface appears to have many holes in it. In fact, only about half of the surface is plotted, and the rest looks like empty space. I know that the integral/function is defined at these empty points, because I can evaluate it there and get a real number. What can I do to fix this? I've tried different integration methods, but still I have no luck. I will upload the code for reference below. Thanks!

test3.mw

> 

with(Statistics):
> 

with(plots):
> 

with(Student[NumericalAnalysis]):
> 

Digits := 3;

3

 (1)
> 

> 

> 

> 

Dmin := 5;

5

 (2)
> 

dosemean := 25;

25

 (3)
> 

dosevar := 500;

500

 (4)
> 

 

 
> 

> 

> 

d1 := RandomVariable(LogNormal(mu, sigma));

_R14171

 (1.1)
> 

eq1 := Mean(d1) = dosemean;

exp(mu+(1/2)*sigma^2) = 25

 (1.2)
> 

eq2 := Variance(d1) = dosevar

exp(2*mu+sigma^2)*(exp(sigma^2)-1) = 500

 (1.3)
> 

distvar := fsolve({eq1, eq2}, {mu, sigma});

{mu = 2.92, sigma = -.767}

 (1.4)
> 

mudist := rhs(distvar[1]);

2.92

 (1.5)
> 

m2 := evalf((1/1000)*round(1000*mudist));

2.92

 (1.6)
> 

``
> 

sigmadist := abs(rhs(distvar[2]));

.767

 (1.7)
> 

s2 := evalf((1/1000)*round(1000*sigmadist));

.767

 (1.8)
> 

``
> 

X := RandomVariable(LogNormal(m2, s2));

_R14172

 (1.9)
> 

Dist := PDF(X, u-Dmin)

piecewise(u < 5, 0, .650*2^(1/2)*exp(-.850*(ln(u-5)-2.92)^2)/((u-5)*Pi^(1/2)))

 (1.10)
> 

``
> 

NN2 := proc (beta, D10f, D10p) options operator, arrow; (1-beta)*exp(-ln(10)*u/D10f)+beta*exp(-ln(10)*u/D10p) end proc;

proc (beta, D10f, D10p) options operator, arrow; (1-beta)*exp(-ln(10)*u/D10f)+beta*exp(-ln(10)*u/D10p) end proc

 (1.11)
> 

NNo := proc (beta, D10f, D10p) options operator, arrow; evalf(Int(Dist*NN2(beta, D10f, D10p), u = Dmin .. 1000, digits = 4, method = _d01ajc)) end proc

proc (beta, D10f, D10p) options operator, arrow; evalf(Int(Dist*NN2(beta, D10f, D10p), u = Dmin .. 1000, digits = 4, method = _d01ajc)) end proc

 (1.12)
> 

eq3 := proc (beta, D10f, D10p) options operator, arrow; NNo(beta, D10f, D10p) = (1-beta)*exp(-ln(10)*REDt/D10f)+beta*exp(-ln(10)*REDt/D10p) end proc;

proc (beta, D10f, D10p) options operator, arrow; NNo(beta, D10f, D10p) = (1-beta)*exp(-ln(10)*REDt/D10f)+beta*exp(-ln(10)*REDt/D10p) end proc

 (1.13)
> 

RED := proc (beta, D10f, D10p) options operator, arrow; fsolve(eq3(beta, D10f, D10p, complex), REDt) end proc;

proc (beta, D10f, D10p) options operator, arrow; fsolve(eq3(beta, D10f, D10p, complex), REDt) end proc

 (1.14)
> 

keq := proc (beta, D10f, D10p) options operator, arrow; ln(1/NNo(beta, D10f, D10p))/RED(beta, D10f, D10p) end proc;

proc (beta, D10f, D10p) options operator, arrow; ln(1/NNo(beta, D10f, D10p))/RED(beta, D10f, D10p) end proc

 (5)
> 

NNo(0.1e-1, 5, 50)

0.4838e-2

 (6)
> 

NormDist := proc (z) options operator, arrow; PDF(RandomVariable(Normal(0, 1)), z) end proc;

proc (z) options operator, arrow; Statistics:-PDF(Statistics:-RandomVariable(Normal(0, 1)), z) end proc

 (7)
> 

``
> 

D1 := proc (z1) options operator, arrow; Dmin+exp(m2+s2*z1) end proc;

proc (z1) options operator, arrow; Dmin+exp(m2+s2*z1) end proc

 (8)
> 

D2 := proc (z1, z2, rho) options operator, arrow; 2*Dmin+exp(m2+s2*z1)+exp(m2+s2*(rho*z1+sqrt(1-rho^2)*z2)) end proc

proc (z1, z2, rho) options operator, arrow; 2*Dmin+exp(m2+s2*z1)+exp(m2+s2*(rho*z1+sqrt(1-rho^2)*z2)) end proc

 (9)
> 

model1 := proc (z1, beta, D10f, D10p) options operator, arrow; (1-beta)*exp(-ln(10)*D1(z1)/D10f)+beta*exp(-ln(10)*D1(z1)/D10p) end proc

proc (z1, beta, D10f, D10p) options operator, arrow; (1-beta)*exp(-ln(10)*D1(z1)/D10f)+beta*exp(-ln(10)*D1(z1)/D10p) end proc

 (10)
> 

model2 := proc (z1, z2, rho, beta, D10f, D10p) options operator, arrow; (1-beta)*exp(-ln(10)*D2(z1, z2, rho)/D10f)+beta*exp(-ln(10)*D2(z1, z2, rho)/D10p) end proc

proc (z1, z2, rho, beta, D10f, D10p) options operator, arrow; (1-beta)*exp(-ln(10)*D2(z1, z2, rho)/D10f)+beta*exp(-ln(10)*D2(z1, z2, rho)/D10p) end proc

 (11)
> 

func1 := proc (z1, beta, D10f, D10p) options operator, arrow; model1(z1, beta, D10f, D10p)*NormDist(z1) end proc;

proc (z1, beta, D10f, D10p) options operator, arrow; model1(z1, beta, D10f, D10p)*NormDist(z1) end proc

 (12)
> 

func2 := proc (z1, z2, rho, beta, D10f, D10p) options operator, arrow; model2(z1, z2, rho, beta, D10f, D10p)*NormDist(z1)*NormDist(z2) end proc;

proc (z1, z2, rho, beta, D10f, D10p) options operator, arrow; model2(z1, z2, rho, beta, D10f, D10p)*NormDist(z1)*NormDist(z2) end proc

 (13)
> 

NNo1 := proc (beta, D10f, D10p) options operator, arrow; evalf(Int(func1(z1, beta, D10f, D10p), z1 = -50 .. 50, method = _d01ajc)) end proc;

proc (beta, D10f, D10p) options operator, arrow; evalf(Int(func1(z1, beta, D10f, D10p), z1 = -50 .. 50, method = _d01ajc)) end proc

 (14)
> 

NNo2 := proc (rho, beta, D10f, D10p) options operator, arrow; evalf(Int(Int(func2(z1, z2, rho, beta, D10f, D10p), z1 = -50 .. 50, method = _d01ajc), z2 = -50 .. 50, method = _d01ajc)) end proc;

proc (rho, beta, D10f, D10p) options operator, arrow; evalf(Int(Int(func2(z1, z2, rho, beta, D10f, D10p), z1 = -50 .. 50, method = _d01ajc), z2 = -50 .. 50, method = _d01ajc)) end proc

 (15)
> 

RED1 := proc (beta, D10f, D10p) options operator, arrow; Dmin-ln(NNo1(beta, D10f, D10p))/keq(beta, D10f, D10p) end proc;

proc (beta, D10f, D10p) options operator, arrow; Dmin-ln(NNo1(beta, D10f, D10p))/keq(beta, D10f, D10p) end proc

 (16)
> 

RED2 := proc (rho, beta, D10f, D10p) options operator, arrow; 2*Dmin-ln(NNo2(rho, beta, D10f, D10p))/keq(beta, D10f, D10p) end proc;

proc (rho, beta, D10f, D10p) options operator, arrow; 2*Dmin-ln(NNo2(rho, beta, D10f, D10p))/keq(beta, D10f, D10p) end proc

 (17)
> 

RAF := proc (rho, beta, D10f, D10p) options operator, arrow; (1/2)*RED2(rho, beta, D10f, D10p)/RED1(beta, D10f, D10p) end proc;

proc (rho, beta, D10f, D10p) options operator, arrow; (1/2)*RED2(rho, beta, D10f, D10p)/RED1(beta, D10f, D10p) end proc

 (18)
> 

evalf(RAF(0, .1, 5, 55))

.717

 (19)
> 

plot('RAF(0, j, 5, 50)', j = 0 .. 1)

 
> 

plot3d('RAF(0, 0.1e-1, p1, p2)', p1 = 1 .. 5, p2 = 10 .. 50);

 
> 

``
> 

NULL
> 

``



Download test3.mw
Please Wait...