Yes, I meant upper case Diff, the inert form

In words to students it would be the same, just in more familiar notation.
And may motivate the more compact form as a variation and explanation.

I do not want to 'hijack' Carl's answer, but:

In my understanding a system parses a command before execution, so this would
be possible. At some costs, of course - but modular arithmetics would need it,
as there are quite often large numbers involved.

Certainly an extra command/function like &^ or pow may be more easy.


But what I miss: that command seems not to be aware of Fermat's little theorem
and an intelligent routine should be able to use it.

While I can accept that modp(f(x), 5) returns f(x) - which now has to be seen
in Z/(5) - I certainly expect that modp(x &^1001, 5) does NOT execute modp at
(`&^`(x,1001), 5) and then return x &^1001 (using the trace command).

I would expect to get x [ to be seen in Z/(5) ] for that case:

x^m = x^n modulo(q) if m = n modulo(q -1) and 1001 = 250*4 + 1 = 1 modulo(5-1),
see http://en.wikipedia.org/wiki/Fermat%27s_little_theorem#Generalizations


Perhaps it should even cover Euler's generalization:

x^m = x^n modulo(k) if m = n modulo( phi(k) ), phi = Euler phi, x coprime to k,
see http://en.wikipedia.org/wiki/Euler%27s_theorem


Though I was not aware of that command (thx, Carl!) it should be improved, IMHO.

That is a poor excuse by Maple. They should implement it ... how to 'overload' it by default?

@acer 

A minor question: you use 64 bit?

It seems to be integer[8] = 8 bytes = "20 places", while I have integer[4] on Maple 32 Bit (and Win 64 Bit) = 4 bytes.

Yes it works - in sense of producing outputs.

But due to the extreme values I would not trust that brute approach (I would go the route already sketched, at least one can not simply feed to the original function - even if I did)

And it ignores the graphics shown by Carl.

I think I have said enough towards that task(s).

@Dima 

I can read and display the file. But it does not execute. And if I copy the text
then it will not execute as well. Do you have a proper sheet as well?

Do you want to say that 

1) there is no k33 because you used k11 + k22 = - k33

2) one has to use S1:=0, S3:=0 and S2:=200 

    ?

  restart:
  approxErf := x -> -Pi^(1/2)/((Pi-1)*x+(x^2+Pi)^(1/2))*exp(-x^2)+1;
  
  L := convert(
       [[1.0, 0.75e-1], [2.0, .1], [3.0, .12], [4.0, .14], [5.0, .16],
       [6.0, .18], [7.0, .2], [8.0, .22], [9.0, .24], [10.0, .26]],
       rational
  ):
  
  NN := -N+(N0B-NB)*(erf((1/2)*x/(sqrt(t)*sqrt(d)))
       -sqrt(erf((1/2)*alpha/d)))/sqrt(erfc((1/2)*alpha/d))+NB:
  
  convert(%, erf):
  aNN:= eval(%, erf=approxErf):
  aNN:= simplify(aNN, power) assuming 0<x, 0<d:
  aNN:= simplify(aNN, size):
  
  alpha:= 2*10^(-12):  NB:= 75/1000:  N0B:= 2/10:  t:= 360000:
   
  NNlog:= eval(NN, d= 10^d):
  aNNlog:= eval(aNN, d= 10^d):
  
  Digits:= 50:
  
  for k to nops(L) do
    af:= numer(simplify(eval(aNN, [x,N]=~ L[k])));
    af_log:= numer(simplify(eval(aNNlog, [x,N]=~ L[k])));
    af_log:= simplify(af_log):
    printf("\nTrying x= %f, N= %f\n", evalf[5](L[k])[]);
    a_dd:= fsolve(af_log, d= -infinity..-0); #10^(-6) .. 10^(-3)
    if a_dd::numeric then
       printf("Solution found: d= %a (for approximation)\n", evalf[15](10^a_dd));
       evalf( eval(af, d=(10^a_dd)) );
       printf("       approx(d) = %a\n", % );
  
       # try to use the 'actual' numerator
       f:= numer(simplify(eval(NN, [x,N]=~ L[k])));
       f_log:= numer(simplify(eval(NNlog, [x,N]=~ L[k])));
       dd:= fsolve(f_log, d= a_dd); #10^(-6) .. 10^(-3)
       if dd::numeric then
         printf("Solution found: d= %a (actual numerator)\n", evalf[15](10^dd));
         evalf( eval(f, d=(10^dd)) );
         printf("            f(d) = %a\n", % );
       else
         printf("No improved solution found.\n");
       end if;
       #
  
    else
      printf("No solution found.\n");
    end if;
  #     print(plots:-semilogplot(af, d= 1e-14..1e-2, labels= [d,'NN']))
  end do:

Trying x= 1.000000, N= 0.075000
Solution found: d= .100000000000000e-29 (for approximation)
       approx(d) = 0.
Solution found: d= .100000000000000e-29 (actual numerator)
            f(d) = 0.

Trying x= 2.000000, N= 0.100000
Solution found: d= .931101154735267e-4 (for approximation)
       approx(d) = .6621e-45
Solution found: d= .864547473017420e-4 (actual numerator)
            f(d) = -.9755e-49

Trying x= 3.000000, N= 0.120000
Solution found: d= .596939616461537e-4 (for approximation)
       approx(d) = .1461e-45
Solution found: d= .570969017295855e-4 (actual numerator)
            f(d) = -.7540e-48

Trying x= 4.000000, N= 0.140000
Solution found: d= .456695783646335e-4 (for approximation)
       approx(d) = .2201e-45
Solution found: d= .445134288859349e-4 (actual numerator)
            f(d) = -.800e-49

Trying x= 5.000000, N= 0.160000
Solution found: d= .355642809533058e-4 (for approximation)
       approx(d) = .2093e-45
Solution found: d= .350843030752839e-4 (actual numerator)
            f(d) = .7242e-48

Trying x= 6.000000, N= 0.180000
Solution found: d= .403074181379548e-10 (for approximation)
       approx(d) = -.1e-48
Solution found: d= .378645553252026e-10 (actual numerator)
            f(d) = .8e-48

Trying x= 7.000000, N= 0.200000
Solution found: d= .100000000000000e-29 (for approximation)
       approx(d) = 0.
Solution found: d= .100000000000000e-29 (actual numerator)
            f(d) = 0.

Trying x= 8.000000, N= 0.220000
Solution found: d= .162872981867713e-18 (for approximation)
       approx(d) = 0.
Solution found: d= .162872981867713e-18 (actual numerator)
            f(d) = 0.

Trying x= 9.000000, N= 0.240000
Solution found: d= .100000000000000e-29 (for approximation)
       approx(d) = 0.
Solution found: d= .100000000000000e-29 (actual numerator)
            f(d) = 0.

Trying x= 10.000000, N= 0.260000
Solution found: d= .100000000000000e-29 (for approximation)
       approx(d) = 0.
Solution found: d= .100000000000000e-29 (actual numerator)
            f(d) = 0.

@Carl Love 

short reply (it is late for me): for the f (above, 27 Jan 2014) I get d ~ .144000000000207e-17

Should I fork to a new answer/question/post or what ever to sketch a way?

Edited to continue:

The idea is as follows (using translation to H:= (x, r) -> erf(x)^2 - erf(x^2/r)
as above, thus looking for large x):

It is easy to see there is a sign change looking at H(r,r), hence there is zero.

For large x that form is messy, because erf quickly equals 1 for numerical input.

So switch to erfc(x) = 1 - erf(x). Now expanding cancels to two 1s and erfc(x)^2
can be (almost) ignored in view of the other terms (care for it later in case of
doubts) and that one expects x to be quite large.

erfc has a simple asymtotics and from that one can expect, that multiplying by
exp(x^2) gives a reasonable scaling (instead of using log). And that does not
change the zeroes.

For moderate r that allows to find a zero, using x0=r as an initial guess.

For general r then I remembered, that for erfc(x)*exp(x^2) ~ Mills ratio there
are inequalities for below and above, which are 'easy' to handle. 

And it is not too hard to find appropriate factors to use that for the expanded
H(x), either using all 3 terms or only 2 (canceling erfc(x)^2).

PS: I would guess that there is only 1 zero.

I find it a bit strange that the OP does not answer, so I keep it short:

It is easy to see that H(r,r) is negative, hence a solution exits. Numerically
it is difficult, since r is very large w.r.t. erf, even if re-written as 'erfc'.

One way might be to look at variations of Mills' ratio and its approximations
by algebraic functions (but I only know those for good absolute errors).

Thus a downward log bisection may be need, though erfc will evaluate to zero.

I think we have r = 6250000/9 with x ~ 694444.444443945 should do it.

Then re-scaling x to the original task has to be done (1st data pair only
others go the same route).

@Carl Love 

So for the first pair we have (rationalized) 

f := d -> 1/8*(-erf(1/1000000000000*1/d)^(1/2)+erf(1/1200*1/(d^(1/2))))/
                erfc(1/1000000000000*1/d)^(1/2)

limit(f(d), d=0, right);
                              -infinity

So together with your experimental plots for RootFinding:-NextZero that says,
there is a zero (if the limit is correct). Which does not say, that one can
find it, I think it is an asymtotic problem (and the poster should say from
where the questions arise).

Rewrite the task: 8*f(d/1440000) becomes g(d, (3/2500)^2) where I use

g := (d, r) -> (-erf(1/(d*r))^(1/2)+erf(1/(d^(1/2))))/erfc(1/(d*r))^(1/2);

We are in case 1 < r. Plot the numerator for r=3/2 and a bit larger too see
what happens and thus reformulate:

h(x,r) = g(1/x^2,r) = erf(x) - sqrt( erf(x^2/r) ) has a zero towards 0, 1<r.

Or more conveniently take H := (x, r) -> erf(x)^2 - erf(x^2/r), since we are
looking for regions with erf(x)^2 < erf(x^2/r). Or other monotonous transforms
staring from that.

While for r=3/2 it is easy (see plot below), it begins to fail numerically beyond
r ~ 7 (as you said due to numerical limitations). 

So I think the graphics are a kind of 'fake' for the limiting situations.


plot for H(x, 3/2):

Edited: hm .., something like erf(x)^2 < erf(x) because this is below 1
and squaring there and erf(x) < erf(x^2/r) for x large by monotony.

@Dima 

"In investment you could find program that I use" - what do you mean by that?

@J4James 

What is your question?

You have model = ... (see above), depending on N, NN, d, x
And you have data L = ...

Now take one data point, say = [1.0, 0.75e-1].

1) What variable of the do you replace by that pair ?
2) What equation do you want to solve for that case ?

Try to type in

convert(A, list);
convert(A, listlist);

The first shows the variables the second with substitutions.

I get the identity matrix, both with Maple 16.02 and 17 concurrent

As expected: Maple substitutes into A.

(Win 7 Pro 64 Bit concurrent , Maple 32 Bit, Intel)

View at the task (cleaned up) in a different way: the Heaviside function is just
an indicator function for a domain. Now try to reduce all to integrate over the
boundary of the domain. Which generally is Stoke's theorem, but here it turns out
to allow a simple case of Green's theorem in the plane.

Files attached:

Part2.mws
Part2.pdf
And in case of the bug in that board:

http://axelvogt.de/temp/Part2.mws
http://axelvogt.de/temp/Part2.pdf