tomleslie

13876 Reputation

20 Badges

15 years, 181 days

MaplePrimes Activity


These are replies submitted by tomleslie

Don't have the FourierSeries package loaded, so can't do a complete check, but you may want to contemplate the limits on integrals in the following


 

f:=t->piecewise(0<=t and t<2*Pi,t^2);

f := proc (t) options operator, arrow; piecewise(0 <= t and t < 2*Pi, t^2) end proc

(1)

a[0]:=int(f(t), t=-Pi..Pi);
a[0]:=int(f(t), t=0..2*Pi);

(1/3)*Pi^3

 

(8/3)*Pi^3

(2)

a[n]:=eval((1/Pi)*int(f(t)*cos(n*t), t=-Pi..Pi), n=2);
a[n]:=eval((1/Pi)*int(f(t)*cos(n*t), t= 0..2*Pi), n=2);

1/2

 

1

(3)

 


 

Download fSeries.mw

answer is infinite - which can be seen from

#
# Set up 'Sum'
#
  mySum:=n^(3/2)*Sum(2*k/(2*k+3), k=1..n);
#
# determine limit of sum as n->infinity
#
  limit(mySum, n=infinity);
#
# Plot the value of sum as n gets very
# large - definitely heading for infinity!
#
  plot(mySum, n=0..1000);

Look at the two terms in your expression, ie n^(3/2) and Sum(2*k/(2*k+3), k=1..n). It is fairly obvious that the first of these terms is headed for infinity as n->infinity. So the only way to get a 'finite' result is if the second terms heads for zero, 'faster' than the first terms for infinity. But look at the terms in this summation: for large 'k' , say k=1000, the value of this term is 2000/2003. In fact for  large 'k', individual terms, asymptote to 1. So as n->infinity, the summation is essentially adding an infinite number of terms, each of which is equal to one. So this summation on its own, would return 'inifinity' as n->infinity.

In other words, both terms in your expresssion -> infinity as n->infinity, and infinity*infinity=infinity

of the distinction betwen local and global optimizers (and generallly don't put too much faiht in the latter!). My point was (and I quote)

In the case of fitting to a sinusoid it really doesn't matter whether you are using 'local' solvers or 'global' solvers. There are always multiple equivalent solutions

I guess those occasions where I recall seeing 'time' measurements expressed in uS (or even nS), I have been using CodeTools:-Usage() with the 'iterations' option.

Using the time() command, I'm tempted to agree that fine resoltuoin is not possible

However with CodeTools:Usage(), I seem to get finer resolution on tiimings.

time() reports the total cputime since the start of the Maple session - rather obviously, this could be a very big number, so what resolution could one reasonably expect?. The sequence

t1:=time()
# some arbitrary maple comands
t2:=time()-t1

is essentially computing the difference between two (not very accurate) numbers. So if one is trying to accurately obtain computation times, then CodeTools:Usage(~someMapleCommand), would seem to be a better alternative - I've had it report timings in microseconds

I posted a longer reply to a similar(ish) comment by 'Acer' on this thread, so I'm only going to give the short version here

Consider the 'phase' variable: if I shift it by a 'period' then I get an identical result. So which of the infinite number of possible values for values for 'phase' did you expect to get???

Similarly frequency: if some value 'f' is a good fit to your data, then 2*f, 3*f..... will provide identically good fits - which one do you want?? And how is a fittingin program supposed to know that what you (probably?) want is

  1. The minimum (postive?) frequency which works
  2. The phase shift should be less than a 'period', where this period is determined by the reciprocal of value of 'frequency' obtained in (1) above

In the case of fitting to a sinusoid it really doesn't matter whether you are using 'local' solvers or 'global' solvers. There are always multiple equivalent solutions

If the function sin(2*Pi*f*t) for  a given value of 'f'' produces a minimal set of residuals, then the function sin(2*Pi*(2*f)*t) will produce an identical set of of residuals - as will 3*f, 4*f etc. So what number for 'frequency' should be returned? f, 2*f, 3*f,......whatever.

This is not an issue of whether one is using 'local' or 'global' - it is simply an observation that any integer muliple of 'frequency' will give exactly the same set of residuals. So which of these infinite number of identical solutions do you want? To the "calibrated human eyeball", the answer may be obvious - the smallest (positive?) one. But how is a fitting program expected to "know" this - as far as such a program is concerned, all integer multiples are equally valid and the one it finally returns will probably depend on whatever initial "guess" was made.

A more or less identical argument can be applied to the 'phase' variable. Assuming (by some miracel) the one has been able to compute the "desired" solution for the frequency variable 'f', then the model equation sin(2*Pi*f*t+phi) will have multiple, identical solutions for the variable 'phi' - just keep adding integer multiples of the period and the solutions are identical. Again the 'calibrated human eyeball' can step in because this 'obviously(?)' requires a 'phase' variable which is less than the period (or 1/frequency). But how does a fitting prgram "know" this? It doesn't. So again, for this parameter, there are an infinite number of solutions which will return an identical set of residuals.

Now the logic of my original post attempts to address these issues, although I'm not sure hiw obvious the thought process is - so I'm going to give it here it very basic terms

  1. Strip out any offset
  2. Having removed any offset, look for zero-crossings. Multiple instances of zero-crossings will enable you to identify a reasonable value for the 'period' of the function, and hence for the frequency, Thus a constraint on the freuqency variable is known
  3. Having obtained an estimete for the 'period' of the function (from (2) above), one now has a constraint on the 'phase shift' - because anything more than a period doesn't count - right!

Mariusz Iwaniuk produced a perfectly reasonable solution for this problem by applying his 'calibrated human eyeball' and coming up with a reasonable set of initial values which would allow the problem to be solved. I have no problem with this.

My original post made an attempt to codify what Mariusz' "calibrated human eyeball" was doing, and hence generate the same solution without introducing some kind of 'magical' initial point.

As my original  post pointed out. my suggested algorithm is not guaranteed to succeed - but it will work most of the time.  Regrettaby, for particularly paranooid input data samples, one may have to revert to the 'calibrated human eyeball' in order to obtain an appropriate set of initial values for the fitting routines

If I set typesetting=extended in Maple 2016 (where the default is typesetting=standard), then the unames() command works as I would expect.

Hence issue seems to be specific to Maple 2017

See my comment on your previous post

Quote

After perusing more of this site I came to the conclusion that one needs to use ones mastery of mathematics and sort it out for ones self (ok, I can manage that)

Actually you do not need to "use ones mastery of mathematics and sort it out for ones self" - the only thing you need to master, is how to upload code to this site, rather than pictures of code. And if you do not know the difference, you are probably beyond help. Try investigating the behaviour of the big green up-arrow in the toolbar.

Alternatively keep your inane comments to yourself - this site is generallly reserved for grown-ups

This is bad (my fault) but also good, because if I change the integration in polar coordinates to allow for the extra factor of r occurring because dx.dy goes to r.dr.theta, the my original response can be rewritten as (change highlighted in red)

   restart:
#
# Define integral
#
   P := (x, y) -> (1/2)*exp(-(1/2)*(x^2+G*y^2-2*B*x*y)/(-B^2+G))/(Pi*sqrt(-B^2+G));
#
# Change to polars
#
   P2:=(r, theta)->changecoords( P(x,y), [x,y], polar, [r,theta]);
#
# Try for symbolic solution in polars - takes ages
# doesn't seem to get anywhere - so comment out
#
#  int(P2(r,theta), theta=0..2*Pi, r=0..infinity);
#
# Try for numeric solution in polars
#
   Int(r*P2(r,theta), theta=0..2*Pi, r=0..infinity);
   evalf(eval(%,[G=2, B=-1]));
#
# Try the "same" integral in cartesians - OOOPs
# different answer
#
   Int(P(x,y), x=-infinity..infinity, y=-infinity..infinity);
   evalf(eval(%,[G=2, B=-1]));

Now, numeric solution in polar and cartesion corrdinates at least agree - and both return 1.0000000, for the supplied values of G,B.

I'm not certain I agree with your comment about integrating only over the scond and fourth quadrants. The OP's final comment is

((int(x = 0 .. infinity))*(int(y = -infinity .. 0))+(int(x = -infinity .. 0))*(int(y = 0 .. infinity)))*P(x, y)

I accept that this is slightly difficult to interpret, because the int() expressions appear to have no integrand - just limits, However since these limits seem to be

x = 0 .. infinity
y = -infinity .. 0
x = -infinity .. 0
y = 0 .. infinity

It looks to me like all four quadrants,  or am I missing something

 

if you listed what you believe to be the "correct" results for these integrals, with an accessible reference, so that I can check

My "simple-minded" approach for evaluating these integrals gives

  1. int(1/ln(t), t = a .. b, AllSolutions) = 0 when a=-1/2, b=1/2
  2. int(sqrt(t^2-1+I*t), t = a .. b, AllSolutions) = -.7808312922-1.178097245*I when a=-1/2, b=1/2. Maple' identifies' the complex coefficient as 3*Pi*(1/8), but provide nothing "useful" for the first term
  3. p3:=int(arctan(t+2*I), t = a .. b, AllSolutions):= 3.673761077*I when a=-1/2, b=1/2, Maple 'identifies' this solution as 7*exp(1)*(1/5)-3*ln(3)*(1/25)

 

So what are the solutions you expect?

which you get is actuall 'Matrix is singular' - and until this is fixed, I don't believe you are going to get anywhere.

You have the most complicated method I have ever seen of producing boundary conditions. Unfortunately it produces incompatilel boundary conditions, in that the same quatity is defined to have two different values. See the final two or three execution groups in the attached. Until this is fixed you have no chance

odesys.mw

You may also want to consider carefully the other highlighted comments I have made in this file.

I apologise for the late response - MaplePrimes has had a problem notifying users of updates to questions, so I was unaware that this discussion was still 'active'

judging by the number of 'notifications' which have turned up for questions I had 'forgotten'. This coud be a long evening :-(

If I just run a quick check on the odes you supply (not sure why you want to treat it as a dynamical system -yet!), then it gets very unstable, very quickly. So I suspect you may be on a loser with this one. Try

    restart
    ode1a := diff(y1(t), t) = 1.052936200*10^5*y1(t)+70106.19000*y2(t)+35169.00000*y3(t);
    ode2a := diff(y2(t), t) = 70106.19000*y1(t)+71031.61000*y2(t)+35511.00000*y3(t);
    ode3a := diff(y3(t), t) = 35169.00000*y1(t)+35511.00000*y2(t)+36100.00000*y3(t);
    sol:=dsolve( [ode1a, ode2a, ode3a, y1(0)=1e-06, y2(0)=1.0e-06, y3(0)=0],numeric);
    plots:-odeplot( sol, [t, y1(t)], t=0..1e-03);
    plots:-odeplot( sol, [t, y1(t)], t=0..1e-03);
    plots:-odeplot( sol, [t, y3(t)], t=0..1e-03);

 

I tried a few different initial conditions, but the solution "exploded" for all ot the somewhere around t=3.0e-2. So unless you have some set of ics for which the above system is 'numerically stable', I don't think that this will ever succeed.

First 112 113 114 115 116 117 118 Last Page 114 of 207