Axel Vogt

5821 Reputation

20 Badges

20 years, 229 days
Munich, Bavaria, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

Literally taken for me the question also makes no sense and as an
ordered finite set is nothing but the according interval of the
Naturals for me that would sufficient (and use substitutions only).

A different source (as I do not have Knuth) and as C program is in
the beautiful online book Jörg Arndt, "Algorithms for programmers",
Chapter 10 at www.jjj.de/fxt/

Have not checked, whether it is Algorithm L.
Thx for the info, but I do not know anything about physical Chemistry :-)
Still I would say, that *your data can not fit the model*, same reasoning.

However I modified your model using 2 location parameters and 

  G := x -> - c
  -a^2*exp((x-mu)*r^2)/(1+exp((x-mu)*r^2))+b^2*exp(x*r^2)/(1+exp(x*r^2))^2;

                        2               2        2        2
                       a  exp((x - mu) r )      b  exp(x r )
        G := x -> -c - -------------------- + ----------------
                                         2                2  2
                       1 + exp((x - mu) r )   (1 + exp(x r ))

From your data and the limiting behaviour at +-infinity of the original
model it follows, that L must be positive and A is negative:

  '[limit(EG(x), x=-infinity),limit(EG(x), x=infinity)]': %: % assuming L<0: 
    '%%%'=%, L<0;
  '[limit(EG(x), x=-infinity),limit(EG(x), x=infinity)]': %: % assuming 0<L: 
    '%%%'=%, 0 < L;

   [      lim         EG(x),      lim       EG(x)] = [A, 0], L < 0
    x -> (-infinity)         x -> infinity


   [      lim         EG(x),      lim       EG(x)] = [0, A], 0 < L
    x -> (-infinity)         x -> infinity

where EG is your original model function, EG:= x -> ((A*exp((2*Pi*x)/L))/
(1+exp((2*Pi*x)/L)))+((B*exp((2*Pi*x)/L))/(1+exp((2*Pi*x)/L))^2);

Trying to fit that always shows a shift off zero: the maximum can not be
achieved close to it, because of A being not 0, since it is responsible
for the tail as well - at that is below -0.02 and certainly lower. 

So I introduced a location parameter mu: some shift along the x-axes for
the "A-term".

Additionally your data at the left side of zero stop at the x-axes.
But I guess that's only true for the provided data and you simply do
not have more. So it seems natural for me to assume data may become
negative as well there. But the limiting behaviour wants them to be
strictly positive.

So I use a level shift parameter c, shifting vertically.

And I came up with the definition for G (close as possible to your model
and the given data):

  EG(x); subs(A=-a^2, B=b^2, L=Pi*2/r^2, %);
  subs(x=x-mu*1,op(1,%)) + subs(x=x-nu*0,op(2,%));
  -c + %;
  G:=unapply(%,x);

Then I get a reasonable fit:

  residuals := map(p -> (G(p[1])-p[2]), data):
  R:= LSSolve(residuals);
  params:=R[2];
  Gfit := unapply(eval(G(x), params),x);
  plot(gfit(x), x=-6 .. 6);
  display([pointplot(data), plot(Gfit(x), x=data[1,1]..data[-1,1])]);

Certainly not superb, but reasonable.

 

will try tonight ... if data go negative to the right that should bear info on signs of A and L, see above

ok, may be ... but if you say the model is correct there still may be the problem of data

Can you at least say what you expect from theory for the wings for large +-x? This never
can be read off the provided data, but all the constallations may be used for a fit, while
only part of them would make sense (i.e. there is no reason, why a minimum through
fitting is a global one and that even exists a global and unique one).

Additionally the model function takes its maximum in x_max := ln(-(A+B)/(A-B))*L/Pi/2
(differentiate f). Your data suggest x_max = epsilon ~ 0 . 

Since L is not zero you will come up with A ~ B/L * Pi * epsilon + ... , epsilon ~ 0. 

That may mean that if you want a close fit over 0 that A is almost 0.

Are you sure for your model function?
I played with that and think there are some issues with your data
(besides certainly the 1st data point should not have a y value of 0).

Let's call f the model function.

Then the limit for large x is 0 respectively A for x= +- infinity and that
depends *only* on the sign of L.

The data however do not allow any guess for that (just looking at the plot),
one does not see the tails (you want to see both of them in best cases).

So you would need enough data to decide that (try the suggested solutions
by substituting L=L^2 or L=-L^2 for new model functions).

The next is f(0) = (A + B/2) / 2, which is responsible for the fit at zero.
This is a bit ugly: it means, that A is given through the tails (by sign L),
while B is given in 0.

Now I was playing more:

  f(x); 'eval(%, B = .1789848-2*A + epsilon*0)'; # 0-value + measurement error
  subs(A=-A^2,%);                                # decide for sign A 
  subs(L=-L^2,%);                                # decide for sign L
  r:=unapply(%,x);
  residuals := map(p -> (r(p[1])-p[2]), data):
  RR:= LSSolve(residuals);
  rfit := eval(r(x), RR[2]);
  plot(rfit(x), x=-6..6);
  display([pointplot(data), plot(rfit(x), x=data[1,1]..data[-1,1])]);

This shows me, that beyond the known data (i.e. the tails) there may be
some which are negative - at least one should know something about that
(i.e. sign A).

For an actual fit I would try B = (2*f_0 - A) + eps (and estimate eps)
and replace 2*Pi*x/L = mu*x (to avoid troubles, even if the posted values
suggest L ~ Pi*sqrt(2)).

If you know your data are strictly positive it is not bad to estimate
log(f(x)) for (x, log(y)) first and feed the result into a second pass
for the actual model.

For a final fit one might think about a weighting to get desired
behaving for peak and tail.

But may be practioneers do know a better way ...

 

Have not followed the whole thread as well, but I think that Maple can do it
(in reply to Thomas Richard and Jacques).

It is more that Maple does not seem to know that 1/GAMMA is holomorphic in the
plane and how to handle 2F1(a,b , m , z) / Gamma(m) for negative integers m.

To avoid typing formulae one often can work with limits and here it is for the
initial data (in Maple 12):

  Int( t^(a-1)*(1-t)^(b-1), t=0..x); value(%);
  subs( x=2, a=3, b=b, % ); value(%); simplify(%);
  limit(%,b=4, right);

which gives -4/5.


Edited to add: have not cared for the usual transformations to the answer given
by Maple (and Euler's integral), but the formula needed is (F = 2F1) 

Limit(F(a,b,c,z)/GAMMA(c),c = -m) = 
  GAMMA(a+m+1)/GAMMA(a)*GAMMA(b+m+1)/GAMMA(b)/GAMMA(m+2) * z^(m+1) * F(a+m+1,b+m+1,m+2,z)  

 

Stumbling is quite natural ... Acer said it quite (anyway I am aware
of my ocassional unpatientience ...)

The question is not silly. To some extend it is ok get a feeling for
the CAS by hands-on and ignoring the example sheets, help or guides.
I understand that quite well, I hated lectures a bit :-)

But from own experience: at least try to follow some systematics, you
will have a better feeling by looking at you learning curve.

So either look what comes with Maple (no, just do it, look at what you
already have and what sounds to interest you) and least work through
the Introductory Programming Guide (ignoring stuff on GUI first) and
the table of contents of Advanced Programming Guide (to have in mind
wher to look things up), both are available for download as pdf.


Besides Maple's application center there are some good (but may older)
lectures on the web. Due to copy rights I do not want to load up such
stuff, but Google may help you, here are some (free for download):


Kraft + "Maple For Math Majors"
Clark + "Symbolic Computations in Mathematics"
Kolinsky + "Einführung in das Computeralgebrasystem Maple"
Geddes + "Symbolic and Numeric Scientific Computation in Maple"
solve(exp(-x)-x = 0);
evalf(%);
identify(%);
                             LambertW(1)
                             0.5671432904
                             LambertW(1)

Most likely somebody may guide you (even if writing in normal style ... stumbling into new is always a bit of pain, do not give up) if you upload your problem

The best I got was Int( P3(xi,m)/2, xi=0 .. x) = 2/5*(V1 - V2) for m < x where

V2 =((-x+2*m)*(x^2-1)^(1/2)+5/2+(x-m)^2)*(x-m)^(1/2)

V1 = (  -I*( -3+m)*EllipticF((x-m)^(1/2)/(m+1)^(1/2)*I,((m+1)/(m-1))^(1/2))
        -I*(3+m^2)*EllipticE((x-m)^(1/2)/(m+1)^(1/2)*I,((m+1)/(m-1))^(1/2)) )*(m-1)^(1/2)

   =  Int(((3/2+1/2*m^2)*x-2*m)/((x-m)*(x^2-1))^(1/2),x)

with V1 and V2 both approaching infinity for x towards infinity

and V1 / V2 = 1 for x ---> infinity (using MultiSeries:-limit)

But then I had to gave up ...
Say your data look like this (with non-constant spacing, but 3 for each line):
 1   15.0000   1.120000E+02
2   15.0003    1.460000E+02
3 15.0006 1.140000E+02
4 15.0008 1.340000E+02
5 15.0011 6.600000E+01
6 15.0014 1.100000E+02
7 15.0017 1.370000E+02
8 15.0019 6.000000E+01
9   15.0022          9.200000E+01
10 15.0025 5.200000E+01
53197 39.9963 9.000000E+01
53198 39.9971 9.900000E+01
53199 39.9978 3.000000E+01
53200 39.9985 4.900000E+01
53201 39.9993 1.330000E+02
and that they are stored in a file data.txt.
Then you write
  fd := fopen(dataFile,READ,TEXT):
  theData:=readdata(fd, 3): # 3 columns
  fclose(fd):
These are 15 tuples:
nData:=nops(theData);
                             nData := 15
Showing them you can use 'print(theData)' or convert(theData, listlist), working with them
may be best by using 'convert(theData, Array)'  or 'convert(%,Array, (datatype = float[8]))'

i refuse to work with those user-unfriendly MW sheets and as they do not even convert to reasonable format (or let allow copy+paste) i just did  C1 + C2 for all parameters  set to one.

which gives 1.

If I get you right, then use variables r1:= ... until your final result and end the lines by a ':' and not by a comma ';'
this prevents Maple from showing the (intermediate) results.
The other way is: just mark the portion of the 'disturbing' sheet using the mouse and through the menu say
Edit/Remove output/From selection

but played with that ... if you disturb a bit P1(0)= c = P2(0) with c = 10^5 - 1 you run into complex solutions and using a series method one sees that for c=10^5 one divides by zero (and have something complex):

'sys'; convert([%],rational): op(%):
G:=dsolve({%, x(0)=1,  P1(0)= c, P2(0) = c, D(x)(0)=0}, fcns, series):

First 76 77 78 79 80 81 82 Last Page 78 of 92