acer

various...

Answered: acer 32343

April 14 2016

1 8

I just posted an answer yesterday to a question about doing this for a BVP.

There I used the numapprox package to get a rational polynomial that approximated the integral well and evaluated very quickly within dsolve/numeric. I suspect that this will allow by far the fastest computation that other approaches (especially with repeated use of solving), and is often reasonably easy to set up.

You could use the same approach for your IVP, I believe. Or you could use a black-box known function approach. And since it's an IVP ODE you could also use the procedure option (but that last one is more complicated to set up.

Here is an IVP example that I just made up. (Here it seems easier than it was for a BVP, when not approximating the integral and allowing the numeric integrations to be deferred until the dsolve/numeric solution is run.) Of course the degree of accuracy that you'll be able to obtain for the numeric integral, and the specifics of the speed performance, will depend on your example (which you didn't post so far).

>

restart;

>	plots:-setoptions(gridlines=false);

>	Digits:=15;

>	Eq:= diff(p(r),r,r) + diff(p(r),r) - f(p(r));

diff(diff(p(r), r), r)+diff(p(r), r)-f(p(r))

>	ICS:=D(p)(0.001)=0, p(0.001) = 3;

For this toy example we can compute the integral exactly. This is just to demonstrate whether the other approaches are accurate.

In general we won't necessarily be able to solve the integral exactly (which is the point of having other approaches).

>	fexpl:=proc(p) local t; #if not p::numeric then return 'procname'(args); end if; int( 3sin(pt), t=0..Pi/2 ); end proc: fexpl(p(r)); fexpl(3.2);

-3*(-1+cos((1/2)*Pi*p(r)))/p(r)

>	Eqexpl:=eval(Eq,f=fexpl);

diff(diff(p(r), r), r)+diff(p(r), r)+3*(-1+cos((1/2)*Pi*p(r)))/p(r)

>	solexpl:=dsolve( {Eqexpl,ICS}, numeric, output=listprocedure ):

>	Pexpl:=eval(p(r), solexpl): Pexpl(1.1);

>	CodeTools:-Usage( plot( Pexpl, 0.001 .. 10.0, size=[600,200]) );

memory used=3.27MiB, alloc change=0 bytes, cpu time=28.00ms, real time=29.00ms, gc time=0ns

Now we'll do it with a procedure which evaluates the integral numerically, on-the-fly. The key thing here is that the procedure fimpl returns unevaluated when passed a symbolic argument such as p(r) for unassigned name r. Note that fimpl(p(r)) is an unevaluated function call in Eqimpl.

This is not as fast. We might be able to speed it up a bit if the numeric integral can be done quicker (at each p value), using options of evalf/Int. But that may not be possible in general, so for the purpose of demonstration I'll leave that alternative commented out.

>

fimpl:=proc(p)
     local t;
     if not p::numeric then return 'procname'(args); end if;

     # This example can be done faster with these additional options for
     # the call to evalf(Int...)).
     # But beware that it may not be able to achieve better accuracy than
     # about epsilon=1e-10 over the full range p=3..4 for this example.
     # And that aspect will be problem specific.
     # Also this limits how well fimpl approximates the integral.
     #      digits=15, epsilon=1e-10, method=_d01ajc

     evalf(Int( 3*sin(p*t), t=0..Pi/2 ));
   end proc:

fimpl(p(r));
fimpl(3.2);

>	Eqimpl:=eval(Eq,f=fimpl);

diff(diff(p(r), r), r)+diff(p(r), r)-fimpl(p(r))

>	solimpl:=dsolve( {Eqimpl,ICS}, numeric, output=listprocedure ):

>	Pimpl:=eval(p(r), solimpl): Pimpl(1.1);

>	CodeTools:-Usage( plot( Pimpl, 0.001 .. 10.0, size=[600,200]) );

memory used=503.04MiB, alloc change=0 bytes, cpu time=3.44s, real time=3.20s, gc time=532.00ms

Now we'll try and speed it up, by incurring an upfront cost of approximating the integral with a rational polynomial in p over the target range of p=3.0..4.5.

The same rationale as given in the comments above apply here as well, with regard to speeding up the numeric integration. Except that here it's a little more convenient as we'll be able to tell before the dsolve/numeric runs whether the integral is computing successfully at all values of p in the desired range.

>	# digits=15, epsilon=1e-10, method=_d01ajc expr:=Int( 3sin(pt), t=0..Pi/2 );

Int(3*sin(p*t), t = 0 .. (1/2)*Pi)

The time to produce the approximating rational polynomial is an up-front cost incurred just once. For this example this up-front cost is even smaller if we use those additional optional arguments in the Int call above.

>	expr_approx:=CodeTools:-Usage( numapprox:-chebpade( expr, p=3..4.5, [8,4] ) ): expr_approx:=eval(expr_approx,T=orthopoly[T]);

memory used=82.84MiB, alloc change=0 bytes, cpu time=596.00ms, real time=545.00ms, gc time=112.00ms

>	plot( expr - expr_approx, p=3..4.5, size=[600,200] );

We need to replace the name p by the function call p(r) in the expression expr_approx.

>	fapprox:=subs(p=p(r),expr_approx);

So now we can construct the new differential equation where f(p(r)) is replaced by our approximating rational polynomial ( in p(r) ).

>	Eqapprox:=eval(Eq, f(p(r))=fapprox);

>	solapprox:=dsolve( {Eqapprox,ICS}, numeric, output=listprocedure ):

>	Papprox:=eval(p(r), solapprox): Papprox(1.1);

The performance of this is now very fast. Mostly it's fast because no numeric integration is now being done when the numeric ODE solution runs (simulates). But it helps for speed to that the approximating expression is just a rational polynomial, as this wouldn't hamper full use of evalhf or the compile option to dsolve/numeric (which I don't do here).

The one-off up-front cost of using numapprox may even be subsumed by runtime savings if we call this dsolve/numeric solution many times (for several plots, etc).

>	CodeTools:-Usage( plot( Papprox, 0.001 .. 10.0, size=[600,200]) );

memory used=2.18MiB, alloc change=0 bytes, cpu time=20.00ms, real time=20.00ms, gc time=0ns

We may also have a look at how the various solutions compare. Here again we can see the timing difference.

>	CodeTools:-Usage( plot( Papprox-Pexpl, 0.001 .. 10.0, size=[600,200]) );

memory used=12.83MiB, alloc change=0 bytes, cpu time=204.00ms, real time=190.00ms, gc time=44.00ms

>	CodeTools:-Usage( plot( Pimpl-Pexpl, 0.001 .. 10.0, size=[600,200]) );

memory used=1.94GiB, alloc change=0 bytes, cpu time=13.92s, real time=12.70s, gc time=2.70s

>

Download ds_numapprox.mw

acer

numapprox...

Answered: acer 32343

April 13 2016

1 2

I executed this in 64bit Maple 2015.2 running on an Intel i5-6600K @ 3.5Ghz. The plots look nicer in Maple itself.

I believe that there would also be a methodology using a blackbox "known" function J. There `diff/J` might also need to be constructed to return a product of two unevaluated function calls, due to the chain-rule. One term in the product would be an unevaluated call to a numeric procedure for the derivative of the integral wrt "f3", and the other a call to represent diff(f3(x),x).

By trying such a "black-box" methodology I was able to get an idea of the solution for f3(x), to know the range for which I needed to approximate the integral. Below I used the numapprox package to approximate the integral over f3(x)=-5..5 but the solution has f3(x)=0..2 or so for x=0..1. Using a narrower range of f3(x) , such as 0..2 say, for computing the approximation of the integral may well be possible. The rational polynomial generated using the numapprox package makes a pretty good approximation and the coefficients are "nice" here.

I was also concerned that I needed to prevent the bvp solver going near f3(x)>14 as there is a singularity of the integral near the right of that -- hence my commented-out use of the approxsoln option to dsolve/numeric (which I may have used, at some point...).

I am not a big fan of using ApproximateInt for this kind of thing. Apart from generating an expression that may need higher working precision in order to avoid round-off error (like the enormous non-nice coefficients the OP obtained) the methods it provides may not approximate the integral all that well. The solution from dsolve/numeric might even appear to compute to small abserr in some cases -- but in a bad case that might just mean it's providing a tight solution to a problem with is "far away" from the original.

Someone else might want to check that this all makes sense.

>	restart; plots:-setoptions(gridlines=false);

>	kernelopts(version);

>	g3 := theta^2*(theta-1)^2; beta := 100; chi := 5; kappa := 5; a := 0;

>	Jexpr := simplify( Int( -beta^2g3/((1-g3f3(x))ln(2kappa(1-g3f3(x)))^2) -chig3/(1-g3f3(x))^4, theta = a .. 1 ), size );

Int(10000*theta^2*(-(1/2000)*ln(10-10*theta^2*(theta-1)^2*f3(x))^2+(-1+theta^2*(theta-1)^2*f3(x))^3)*(theta-1)^2/((-1+theta^2*(theta-1)^2*f3(x))^4*ln(10-10*theta^2*(theta-1)^2*f3(x))^2), theta = 0 .. 1)

>

dsys3[3] := -0.104166022802734e-3*(diff(f3(x), x, x, x, x, x, x))
            + 0.272050542507459e-1*(diff(f3(x), x, x, x, x))
            - 2.43720057325002*(diff(f3(x), x, x))-.16558927348305*(diff(f1(x), x, x, x))
            + 1.17598588657475*(diff(f1(x), x))+2.73022914514365*10^(-10)*(diff(f2(x), x, x))
            - 3.36155011672494*10^(-9)*f2(x)+1.98218022588292*f3(x) + Jexpr;

>	Digits:=15: Japprox := CodeTools:-Usage( numapprox:-chebpade(subs(f3(x)=f3x,Jexpr), f3x=-5..5, [6,4]) );

memory used=0.77GiB, alloc change=460.18MiB, cpu time=6.10s, real time=5.83s, gc time=628.00ms

>	Japprox := eval(Japprox, T=orthopoly[T]);

(-53.8423993376853+5.75565686683418*f3x-.177634731674162*f3x^2+0.138034572869404e-2*f3x^3+0.204298775667079e-5*f3x^4+0.911892344321782e-8*f3x^5+0.333247034619404e-10*f3x^6)/(.854136473658405-.167536306557109*f3x+0.116105495497653e-1*f3x^2-0.329366191542557e-3*f3x^3+0.312173640332461e-5*f3x^4)

>	Digits:=15: plot( subs(f3(x)=f3x,Jexpr) - Japprox, f3x=-5..5, size=[600,100], adaptive=false, numpoints=50 );

>

dsys3[1] := 1.39804733500615*(diff(f1(x), x, x, x, x))-20.8373457667114*(diff(f1(x), x, x))
            - 2.82678993247240*(diff(f2(x), x, x, x))+52.7168505310530*(diff(f2(x), x))
            + 0.662357093932194e-1*(diff(f3(x), x, x, x))
            - .470397869798952*(diff(f3(x), x))-3.30648328006403*f1(x):

>

dsys3[2] := 8.81246609032610*(diff(f2(x), x, x, x, x))-365.238395383402*(diff(f2(x), x, x))
            + 935.794579019693*f2(x)+2.82693397272729*(diff(f1(x), x, x, x))
            - 30.1911843713272*(diff(f1(x), x))+1.09281316666667*10^(-10)*(diff(f3(x), x, x))
            - 1.34509666666667*10^(-9)*f3(x):

>	dsys3[3] := subs( Jexpr=subs(f3x=f3(x),Japprox), dsys3[3] );

>

dsys := {dsys3[1], dsys3[2], dsys3[3],
         f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0,
         ((D@@1)(f1))(0) = 0, ((D@@1)(f1))(1) = 0,
         ((D@@1)(f2))(0) = 0, ((D@@1)(f2))(1) = 0,
         ((D@@1)(f3))(0) = 0, ((D@@1)(f3))(1) = 0,
         ((D@@2)(f3))(0) = 0, ((D@@2)(f3))(1) = 0}:

>

Digits:=15:
sol:=CodeTools:-Usage(
       dsolve(dsys, numeric
              #, maxmesh = 512, abserr = 1e-7 ## works with Digits:=18
              , maxmesh = 256, abserr = 1e-5 ## works with Digits:=15
              , method = bvp[middefer]
              , range = 0 .. 1
              , output = listprocedure
              #, approxsoln = [f3(x)=1/2, f1(x)=x, f2(x)=x]
             ) ):

memory used=8.42MiB, alloc change=1.37MiB, cpu time=480.00ms, real time=427.00ms, gc time=100.00ms

>	plot(eval(x,sol), 0..1, size=[600,100],color=black);

>	plot(eval(f1(x),sol), 0..1, size=[600,100]); plot(eval(f2(x),sol), 0..1, size=[600,100]); plot(eval(f3(x),sol), 0..1, size=[600,100]);

>	plot(eval(diff(f1(x),x),sol), 0..1, size=[600,100], color=green); plot(eval(diff(f2(x),x),sol), 0..1, size=[600,100], color=green); plot(eval(diff(f3(x),x),sol), 0..1, size=[600,100], color=green);

>	plot(eval(diff(f3(x),x,x),sol), 0..1, size=[600,100], color=blue);

>

Download bvp_numapprox.mw

acer

set component values only at after the l...

Answered: acer 32343

April 09 2016

1 0

Setting and resetting the component values repeatedly, while the loop runs to find your pair of integers, may not be what you want.

If that's the case then I suggest that you write the loop to find the x and y values (no % prefixes, etc). And then, after the loop finishes, make a single pair of calls to DocumentTools:-Do or DocumentTools:-SetProperty in order to put those values into the two components.

acer

rationals...

Answered: acer 32343

April 06 2016

0 1

The easiest way out of that mess is to not use floats so early on in your computations. Create your expressions using exact rationals instead of floats. If you feel an inescapable desire to assign values to names before creating your expressions (ie, assigning values to parameter names) then assign exact rational values. That's especially important if you want to do symbolic manipulation or simplification of the expressions, and your problem here is not atypical.

On the other hand, if you got handed this problematic expression as the output of some command, then you could try something like one of these (which work for this example but may need tweaking if you have more or fewer trailing digits).

Let's suppose your expression were assigned to the name `ee`.

simplify(identify(evalf[9](ee)));

simplify( convert(evalf[9](ee), rational, exact) );

Note the the former may also resolve things like a float approximation to sqrt(2). And note that the latter produces a result that has no floats (which may be a Good Thing).

The tricky bit is the evalf[9], as the amount of rounding you may need will be problem specific.

If you ever need to work with floats you can also apply the evalf command to your "exact" expression. But hold off until you need to. It's usually far easier to convert from exact to floats accurately than it is to convert back robustly.

acer

precision and accuracy...

Answered: acer 32343

April 06 2016

1 1

The main reason for your problem is that your numeric integral is being computed too coarsely for fsolve to be able to achieve its target accuracy (which depends on the value of Digits).

fsolve will try and obtain approximately 10^(-Digits) degree of accuracy, and it can only obtain that if the expression it is root-solving converges to that many digits. You have Digits set to 15.

You had the evalf/Int call set to obtain only epsilon=1e-5 accuracy. So the thing that fsolve was root-finding was returning something which could vary (effectively randomly) beyond the 5th decimal.

So fsolve would never be able to attain a root accurate to 10(-Digits), because the thing it was root-finding was far less accurate than that.

The solution is to control the accuracy and working precision better. You didn't write how accurate you want the root. But your explicit assignment of Digits:=15 at the top of your sheet affects the accuracy that fsolve is trying for.

In practice the _cuhre method will be able to get about 1e-13 degree of accuracy in the results for your integrand over the given range with some reliability. So you'll need the value of Digits (at the level at which you call fsolve) to be no more than 13. You'll get better performance here if you set the tolerance for the numeric integration to be one digit finer than what fsolve is expecting. And you'll get better results from the numeric integration if you set the working precision for evalf/Int to be about two digits greater than what is used for the (negative of the degree of the) tolerance.

So I'll attempt to get Digits as high as I can, while keeping the evalf/Int tolerance as coarse as I can (since that affects the speed of the whole computation). You can see from the results below that you'll get success much faster if you are willing to accept a root that is less accurate. But even at Digits=10 and ten correct places in the root it is pretty fast (less than two seconds on my machine).

I'm using 64bit Maple 2015.2 on Linux. Results may vary with platform.

>

restart:
Digits:=15:
exprr1:=0.502654861826898e-15*(0.148742399611976343441959840319595e-3*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(3, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.148742399611976343441959840319595e-3*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(3, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.312166710257963624050280526325985e-3*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(2, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(2, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.312166710257963624050280526325985e-3*LaguerreL(2, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(2, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.710593209089338584778833390642454e-3*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(2, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.206721610191267713266886470902853e-3*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(3, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.206721610191267713266886470902853e-3*LaguerreL(3, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.143690143152173589759516081260642e-2*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(2, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.143690143152173589759516081260642e-2*LaguerreL(2, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.398435484547028829404909401385469e-3*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(4, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.398435484547028829404909401385469e-3*LaguerreL(4, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+.107493033305967349851691237358591*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.133119364239157944745420092022615e-1*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.405446986991729344149682899708784e-1*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.405446986991729344149682899708784e-1*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.153908710106461097713933570890594e-2*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(2, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.898474989760274443908140074783779e-2*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.898474989760274443908140074783779e-2*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.249116576897628560729046248654299e-1*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(2, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.249116576897628560729046248654299e-1*LaguerreL(2, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.419228231638236639352049290560905e-1*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.313119716214743723056584711370644e-3*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(3, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.140850157341145270857085560645989e-2*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(2, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.140850157341145270857085560645989e-2*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(2, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.356116251994039684832785446262953e-2*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(2, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.236594125595324402353664739536234e-2*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(3, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.236594125595324402353664739536234e-2*LaguerreL(3, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.283204531922458322873194261111897e-2*LaguerreL(2, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(2, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.356116251994039684832785446262953e-2*LaguerreL(2, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.672871032108499739917557182214943e-2*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(1, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.571268576092047717830863567470608e-2*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(3, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.571268576092047717830863567470608e-2*LaguerreL(3, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.144041670766901641594252874691171e-1*LaguerreL(1, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(2, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))-0.144041670766901641594252874691171e-1*LaguerreL(2, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(1, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(0, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2))+0.557948633363586e-4*LaguerreL(0, .794726420679134*r2+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.58902025541558*r2*t-1.58902025541558*d1val)*LaguerreL(0, 1.58902025541558*r2*t+1.58902025541558*d1val+0.635781183662037e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-.794726420679134*r2)*LaguerreL(4, 1.84574775857457*r2+3.69049083866691*r2*t+3.69049083866691*d1val-0.147659831629265e-15*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)))^2*(exp(-.922873879287283*r2+0.102517974484287e-16*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)-1.84524541933345*r2*t-1.84524541933345*d1val))^2*(1.99945567942447*r2*t+1.99945567942447*d1val)*r2*(0.624999953680278e33*d1val^2-0.850963970396069e29*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.850500773197947e29*r2^2)^(1/2)/(d1val*(.249999981472111*d1val^2-0.340385588158428e-4*(1.99945567942447*r2*t+1.99945567942447*d1val)^2+0.340200309279179e-4*r2^2)^(1/2)):

>	a1:=.49986388281: a2:=.49986388281: a3:=.0002722343: ### This part above is just a very small input

>

>	### This part simplifies the above input expr0:=combine(exprr1,power): expr1:=simplify(%,sqrt): expr1:=simplify(%,LaguerreL): eee:=(expr1*r2)/(a2+a3):

>	### This is the integration step

>	drval:=subs(igrand=eee, proc(d1val,dig,eps) #print(Digits,dig,eps); evalf(Int(Int(igrand,t=-a2..a2),r2=0..d1val/a2, ':-digits'=dig,epsilon=eps,method=_cuhre)); end proc):

>	Logdr:=proc(r, dig, eps)

>	local val: # Increase working precision to avoid round-off in the igrand # evaluations. Think of this as "guard digits" for igrand. # Note that this does not affect `dig` which is passed for the # working precision of evalf/Int's method's control code (cuhre). Digits:=max(dig, Digits+3);

>	val:=drval(r, dig, eps):

>	return val:

>

end proc:

>	### This finds the peak position and the value at this position Optedd:=Optimization[NLPSolve](evaln(Logdr(var,12,1e-5)),var=0..3, optimalitytolerance=3E-15, assume = nonnegative,maximize=true): peak_position:=solve(Optedd[2][1],var); peak_max:=evalf(Optedd[1]);

>	### This calculates the right hand value of the FWHM calculation. Digits:=9; eps:=Digits+1; prec_for_int:=min(15,eps+2); forget(evalf); forget(`evalf/int`); forget(fsolve); FWHMR:=CodeTools:-Usage( fsolve(evaln(Logdr(var,prec_for_int,1.0*10^(-eps)))-peak_max/2, var=peak_position..3) );

memory used=4.24MiB, alloc change=0 bytes, cpu time=1.05s, real time=1.05s, gc time=0ns

>	Digits:=10; eps:=Digits+1; prec_for_int:=min(15,eps+2); forget(evalf); forget(`evalf/int`); forget(fsolve); FWHMR:=CodeTools:-Usage( fsolve(evaln(Logdr(var,prec_for_int,1.0*10^(-eps)))-peak_max/2, var=peak_position..3) );

>

memory used=3.83MiB, alloc change=32.00MiB, cpu time=1.44s, real time=1.44s, gc time=16.00ms

>	Digits:=11; eps:=Digits+1; prec_for_int:=min(15,eps+2); forget(evalf); forget(`evalf/int`); forget(fsolve); FWHMR:=CodeTools:-Usage( fsolve(evaln(Logdr(var,prec_for_int,1.0*10^(-eps)))-peak_max/2, var=peak_position..3) );

memory used=3.98MiB, alloc change=0 bytes, cpu time=1.94s, real time=1.94s, gc time=0ns

>	Digits:=12; eps:=Digits+1; prec_for_int:=min(15,eps+2); forget(evalf); forget(`evalf/int`); forget(fsolve); FWHMR:=CodeTools:-Usage( fsolve(evaln(Logdr(var,prec_for_int,1.0*10^(-eps)))-peak_max/2, var=peak_position..3) );

memory used=3.68MiB, alloc change=2.00MiB, cpu time=2.65s, real time=2.66s, gc time=0ns

>	Digits:=13; eps:=Digits+1; prec_for_int:=min(15,eps+2); forget(evalf); forget(`evalf/int`); forget(fsolve); FWHMR:=CodeTools:-Usage( fsolve(evaln(Logdr(var,prec_for_int,1.0*10^(-eps)))-peak_max/2, var=peak_position..3) );

memory used=12.05MiB, alloc change=0 bytes, cpu time=12.79s, real time=12.81s, gc time=0ns

>

# Here it breaks down, unless we allow eps to be the same degree as Digits. This is
# pushing it, as the results coming out of evalf/Int are going to be no finer than fsolve's
# target accuracy. This can take a long time.
Digits:=14;
eps:=Digits;
prec_for_int:=min(15,eps+2);
forget(evalf); forget(`evalf/int`); forget(fsolve);
FWHMR:=CodeTools:-Usage( fsolve(evaln(Logdr(var,prec_for_int,1.0*10^(-eps)))-peak_max/2,
var=peak_position..3) );

memory used=66.30MiB, alloc change=0 bytes, cpu time=70.46s, real time=70.52s, gc time=24.00ms

>

# Again it breaks down. The maximal allowed setting for prec_for_int is 15,
# and the minimal allowed `epsilon` is 0.5*10^(1-15). But that's not going to
# produce a result accurate enough for fsolve to be satisfied about convergence
# with Digits=15 at this top level. This may never succeed.
Digits:=15;
eps:=14.3; 10^(-eps);
prec_for_int:=round(min(15,eps+2));
forget(evalf); forget(`evalf/int`); forget(fsolve);
FWHMR:=CodeTools:-Usage( fsolve(evaln(Logdr(var,prec_for_int,1.0*10^(-eps)))-peak_max/2,
var=peak_position..3) );

memory used=22.56MiB, alloc change=32.00MiB, cpu time=31.91s, real time=31.92s, gc time=24.00ms

>	# It might be tough to get the root more accurate than this with the cuhre method.

Download fsolve_numl_int_modif.mw

I suggest that you do not use DirectSearch for this problem, unless you don't care about knowing how accurate the root is.

acer

double precision machine constants...

Answered: acer 32343

April 06 2016

0 0

These may vary (a little) amongst the supported operating systems.

evalhf(DBL_MIN);
                                               -307
                         1.00000010000000008 10    
evalhf(DBL_MAX);
                                               307
                          9.9999999000000001 10   
evalhf(DBL_EPSILON);
                                                -16
                          2.22044604925031308 10

acer

symmetric wrt names...

Answered: acer 32343

April 05 2016

0 0

I prefer a scheme which appears "symmetric", ie. not to prefer one name over another, and not to have to patch up any term referenced by position or on a scheme not robust to renaming of all the variables in play, etc.

restart:

ee:=(-(1/2)*ib-(1/2)*ia+ic)*vc+(-(1/2)*ia-(1/2)*ic+ib)*vb+(-(1/2)*ic-(1/2)*ib+ia)*va:

rule:=ia+ib+ic=0:

eqs:=map(-rhs=-lhs,map2(isolate,rule,indets(rule,name)));

                            eqs := {ia + ib = -ic, ia + ic = -ib, ib + ic = -ia}

subs(eqs,map(simplify,ee,{ia+ib+ic=0}));

                                        3 vc ic   3 vb ib   3 va ia
                                        ------- + ------- + -------
                                           2         2         2

Obvious attempts using algsubs or simplify-with-side-relations seem to produce a result which, while equivalent to the desired target, contains four terms. This can be patched up by repairing just the offending pair of terms, it seems like "cheating" a little to me if this has to re-treat two of the terms while holding the others static. If all the names were changed to something else, breaking their lexicographic order, then the simplification scheme should ideally not need to be changed. Just some thoughts.

collect(simplify(ee, {rule}),[va,vb,vc]);      

                                   3 va ia   3 vb ib      /  3 ia   3 ib\
                                   ------- + ------- + vc |- ---- - ----|
                                      2         2         \   2      2  /

acer

multi-line caption...

Answered: acer 32343

April 03 2016

2 0

There are a few differences between how title and caption (for 3D plots) are handled versus how legend (for 2D plots) are handled.

title and caption can't be placed on the left or right
plots:-display uses only the last caption title from a sequence of plot arguments, while it merges legends
title and caption don't automatically get an graphic icon for the object, as legends do
typeset cannot be applied to math expressions for titles and captions, as it can be for legends. But instead Typesetting:-Typeset can be applied

But you can construct a multi-line caption (or title) for a 3D plot which appears like a legend. We can build the graphic icons for the surface/curve, partly by using techniques as mentioned here.

restart;

with(Typesetting): with(plots):

clist := ["Cyan", "Orange"]:

P1 := plot3d(sin(x)*(1-y)^2, x=-10..10, y=-1..1, color=clist[1]):
P2 := plot3d(cos(x)*(1-y)^2, x=-10..10, y=-1..1, color=clist[2]):

lg := mrow(mn("        ",mathbackground=clist[1]),
                     mn("  "), Typesetting:-Typeset(sin(x)*(1-y)^2),
                     mn("\n        ",mathbackground=clist[2]),
                     mn("  "), Typesetting:-Typeset(cos(x)*(1-y)^2)):

display(P1, P2, title=lg);

display(P1, P2, caption=lg);

restart;

with(ColorTools): with(Typesetting): with(plots):

P1 := plot3d(sin(x)*(1-y)^2, x=-10..10, y=-1..1, color="SkyBlue"):
P2 := plot3d(cos(x)*(1-y)^2, x=-10..10, y=-1..1, color="PeachPuff"):
P3 := spacecurve([t,1/2,cos(t)+3], t=-10..10, thickness=2, color="Green"):
P4 := plot3d(sin(x)*(1-y)^2, x=-10..10, y=-1..-1/2, color="Red",
             style=point, grid=[30,10], symbol=solidsphere):

#
# Picking off the color could be done with a procedure, which
# naturally would be careful about the various formats that the
# COLOR substructure might take.
#
select(type, [op([1,..],P1)], 'specfunc(anything,{:-COLOR,:-COLOUR})')[1]:
cP1 := RGB24ToHex(ToRGB24(Color([op([2..],%)]))):
select(type, [op([1,..],P2)], 'specfunc(anything,{:-COLOR,:-COLOUR})')[1]:
cP2 := RGB24ToHex(ToRGB24(Color([op([2..],%)]))):
select(type, [op([1,..],P3)], 'specfunc(anything,{:-COLOR,:-COLOUR})')[1]:
cP3 := RGB24ToHex(ToRGB24(Color([op([2..],%)]))):
select(type, [op([1,..],P4)], 'specfunc(anything,{:-COLOR,:-COLOUR})')[1]:
cP4 := RGB24ToHex(ToRGB24(Color([op([2..],%)]))):
lg := mrow(

           mn("\n           ",mathbackground=cP1, size=10),
           mn("  "), Typesetting:-Typeset(sin(x)*(1-y)^2),

           mn("\n           ",mathbackground=cP2, size=10),
           mn("  "), Typesetting:-Typeset(cos(x)*(1-y)^2),

           mn("\n     "), mn(cat(" "," "$42," "),mathbackground=cP3, size=3),
           mn("  "), Typesetting:-Typeset([t,0.5,cos(t)+3]),

           mn(cat("\n","&CenterDot; "$3),mathcolor=cP4, size=14),
           mn("  "), Typesetting:-Typeset(sin(x)*(1-y)^2)

           ):

plots:-display(P1, P2, P3, P4, caption=lg);

It should even be possible to to write procedures which 1) construct such captions from what is found in a PLOT3D structure (though any math expression would have to be passed separately), and 2) merge captions from separate 3D plots into a multi-line caption, while also merging the plots.

If desired the results from those calls to Typesetting:-Typeset could get size modification, to use a smaller math font. size, I suspect.

acer

surfdata (Maple 2015)...

Answered: acer 32343

April 03 2016

1 4

>

restart:

>	with(plots):

>	nx:=20: tmax:=50: T1:=1: T2:=10: L:=1: k:=1: rho:=1: cp:=1: chi:=k/rho/cp: h:=L/(nx-1): t:=1e-3:

>	T:=Array(0..nx,0..tmax,datatype=float[8]): for k from 0 to nx do T[k,0]:=T1 od:

>	for w from 1 to tmax do

>	T[0,w]:=T1: T[nx,w]:=T2:

>	for q from 1 to nx-1 do

>	T[q,w]:=T[q,w-1]+chit/h^2(T[q+1,w-1]+T[q-1,w-1]-2*T[q,w-1]);

>

od:

>

od:

>

surfdata(Matrix(T), 0..nx, 0..tmax, style=surface,

         dimension=2, labels=["x", "t"],
         #dimension=3, labels=["x", "t", "T"],

         gridsize=[100,100], # interpolate, for smoothness

         colorscheme=["zgradient",
                      ["Indigo","Blue","Cyan","Green","Yellow","Orange","Red"]]);

Download tempdens.mw

acer

parametric infinite sums...

Answered: acer 32343

April 02 2016

1 0

See the help topic updates/Maple2016/compatibility in your Maple 2016.

Download formalsum2016.mw

acer

IterativeMaps...

Answered: acer 32343

March 31 2016

1 3

(It turns out that I mostly just had to look for the file where I'd done pretty much this, earlier.)

If the IterativeMaps commands also give you grief you could also implement it (slower) using plot3d or densityplot.

As mentioned in the code below, increasing `maxiter` will reduce the areas of nonconvergence (black areas) for such univariate polynomials.

I used Maple 2015.2 on 64bit Linux.

>

restart:

>	with(IterativeMaps): with(ImageTools):

>	ee := z^5 - 1;

>	## For ee a polynomial we can compute all the roots. ## This is used to specify the colors, below. rts:=[fsolve(ee,z,complex)]: evalf[4](%);

>	minr,maxr := (min,max)(Re~(rts)): mini,maxi := (min,max)(Im~(rts)):

>	dee := diff(ee, z): Z := eval( z - ee/dee, z=zr+zi*I):

>	fzi := evalc( Im( Z ) ):

>	fzr := evalc( Re( subs(zi=zit, Z) ) ):

>	N := 400; maxiter := 40; # increasing this will reduce nonconvergence h,w := N, floor(N*(maxr-minr)/(maxi-mini));

>

newt := Escape( height=h, width=w,
                       [zi, zr, zit],
                       [fzi, fzr, zi],
                       [y, x, y], evalc(abs(eval(ee,[z=zr+zit*I])))^2<1e-10,
                       2*minr, 2*maxr, 2*mini, 2*maxi,
                       iterations = maxiter,
                       redexpression = zr, greenexpression = zit, blueexpression=i ):

>	## View the iteration count layer #Embed([ FitIntensity(newt[..,..,3]) ]);

>	## By default use full saturation and full intensity img:=Array(1..h,1..w,1..3,(i,j,k)->1.0,datatype=float[8]):

>	## Optional, set intensity to zero for non-convergence img[..,..,3]:=Threshold(newt[..,..,3],1,low=maxiter+1): img[..,..,3]:=Threshold(img[..,..,3],maxiter-1,high=0,low=1):

>	## Optional, scale the saturation by the time to converge img[..,..,2]:=1.0-~FitIntensity(newt[..,..,3]):

>	## Done inefficiently, set the hue by position in the list of ## all roots G:=Array(zip((a,b)->min[index](map[evalhf](abs,rts-~(a+bI))), newt[..,..,1],newt[..,..,2]),datatype=float[8],order=C_order): img[..,..,1]:=240FitIntensity(G):

>	## Convert from HSV to something viewable final := HSVtoRGB(img):

>	Embed(final);

Download newtbasin0.mw

And here it is for N=600 and maxiter=200. Using a logarithm for adjusting the saturation layer might be a way to get a nicer impression of the convergence rate.

acer

op...

Answered: acer 32343

March 30 2016

1 2

c:=INTERVAL(1.41421356167 .. 1.41421356308):

uo,low:=op(op(c));                          

                           uo, low := 1.41421356167, 1.41421356308

Or,

up,low:=op([1,..],c);

                           up, low := 1.41421356167, 1.41421356308

acer

add, if...

Answered: acer 32343

March 29 2016

0 0

You can use the add and the if command together.

Note the single backticks (left-quotes) around the word if.

V := Vector[row]( 10, i -> i^2 );

                 V := [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]

add( `if`( v>30, v, 0), v=V );

                                     330

add( map( v -> `if`( v>30, v, 0 ), V ) );

                                     330

add( piecewise( v>30, v, 0 ), v=V );

                                     330

M := Matrix( 5, 5, (i,j) -> (i+j) );

                                 [2  3  4  5   6]
                                 [              ]
                                 [3  4  5  6   7]
                                 [              ]
                            M := [4  5  6  7   8]
                                 [              ]
                                 [5  6  7  8   9]
                                 [              ]
                                 [6  7  8  9  10]

add( `if`( M[i,3] > 5, M[i,3], 0), i=1..5 );

                                     21

add( `if`( v > 5, v, 0), v = M[..,3] );

                                     21

Is that the kind of thing you mean?

acer

manually (in a few ways)...

Answered: acer 32343

March 28 2016

1 0

>

restart;

>

#
# One way to consider the problem is to notice that the subterm
# 18*s^2+32 happens to equal 9*s^2 + 9*s^2+32 where those
# two parts are the M^2 and N^2 of the desired (M+N)^2 target.
#
# We can forcibly separate those enough to allow factor or CompleteSquare
# to obtain the desired form. But how to automate the split?
#

>	a := 6ssqrt(9s^2 + 32) + 18s^2 + 32;

>	t := indets(a, `^`(anything, 1/2))[1];

>	ft := freeze(t);

>	subs(t = ft, a - t^2 + ft^2); # a rather manual step

>	thaw(factor(%));

>	# But how do we know how to divvy it up? How to find A=1? subs(t = ft, a + A*(-t^2 + ft^2)); Student:-Precalculus:-CompleteSquare(%, ft); subs(S=s^2, collect(algsubs(s^2=S, %), [S], u->simplify(u, size))); thaw( eval(%, A=1) );

>	restart; # Once again, if you know in advance how to split terms, it's easy. # How to automate it? eq1 := 6ssqrt(9s^2 + 32)=2MN: eq2 := 18s^2 + 32=M^2+N^2: map(factor, eq1+eq2); eval(%, solve( {eq1,eq2}, {M,N}, explicit )[1]);

>

Download factorfun.mw

acer

augmenting the context-menu...

Answered: acer 32343

March 25 2016

2 5

If the following code is saved to a Maple initialization file then a new entry for calling the RationalFunctionPlot command should appear in the right context-menu, under the existing "Plots" submenu, if the expression is of type `ratpoly`.

You can of course adjust this. The `caption` might not be set to the empty string. Or several such new menu items could be added, in their own subsubmenu, to provided various common choices for the command's options.

:-__AugmentCM := :-ContextMenu:-CurrentContext:-Copy():
:-__AugmentCM:-Entries:-Add("Rational Function",
            ":-Student:-Precalculus:-RationalFunctionPlot(%EXPR,':-caption'=\"\")",
            ':-ratpoly',
            ':-operator'=:-Typesetting:-mover(Typesetting:-mo("→"),
                         :-Typesetting:-mi("Rational Function Plot")),
            ':-category'="Category 1",
            ':-helpstring'="Rational Function Plot",
            ':-submenu'=["Plots"]
           ):
:-ContextMenu:-Install(:-__AugmentCM);

Or you could adjust this code to invoke instead the Student:-Precalculus:-RationalFunctionTutor command. This command is already available via the main menubar choice Tools->Tutors->Precalculus but invoking it from there won't automatically pass in the desired expression. But if this command is added to the context-menu then the right-clicked expression will be passed in and used.

See the help page with topic worksheet,reference,initialization for documentation on how to create an initialization file.

If you want all the students to get this enhancement automatically for a network installation of Maple then you just add it to the (possibly new) file <MAPLE>/lib/maple.ini where <MAPLE> is the location of the Maple installation. Or do it for each machine, if Maple is installed locally on machines. Test it first, of course.

acer

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

32343 Reputation

29 Badges

Social Networks and Content at Maplesoft.com

MaplePrimes Activity

These are answers submitted by acer

various...

numapprox...

set component values only at after the l...

rationals...

precision and accuracy...

double precision machine constants...

symmetric wrt names...

multi-line caption...

surfdata (Maple 2015)...

parametric infinite sums...

IterativeMaps...

op...

add, if...

manually (in a few ways)...

augmenting the context-menu...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again