tomleslie

13876 Reputation

20 Badges

15 years, 181 days

MaplePrimes Activity


These are replies submitted by tomleslie

You used the term "Rayleigh Identity" (which I'd never heard of) in you original worksheet. When I looked at the sum you were actually performing, it seemed (to me) to be a verification of (what I think of as) Parseval's equation. This was (sort of) confirmed by Wikipedia - so I really didn't give the matter much more thought. The distinction between the two didn't seem to be relevant for your problem.

My understanding is that (strictly speaking) Parseval's equation applies to square-integrable functions and their Fourier transform. (NB, not Fourier series) Since it involves integrals from -Infinity..Infinity in both time and frequency domains, tt can be interpreted as a statement that the total "energy" of the signal measured in the time-domain must be the same as the total "energy" measured in frequency domain.

Note that, in this form, Parseval's identity cannot be applied to periodic functions because in general these are not square-integrable from -Infinity to Infinity - eg the integral of sin(t)^2 from -Infinity to Infinity is infinite.

So I guess that Rayleigh's Identity is a statement about periodic functions and their Fourier series which is (sort of) equivalent to Parseval's theorem for square-integrable functions and their Foutier transform. As a practical matter, by restricting the time-domain integral to a single period, the integration can be performed. My interpretation of this is that Parseval's theorem is a statement about signal "energy", which would be infinite for a periodic function. Restricting the time-domain integral to a single period means that Rayleigh's identity (in some sense) is a statement about signal "power" (ie energy/time) rather than energy.

So I probably shouldn't have used the term "Parseval" anywhere in my response to your post.

My only excuse is that I (like some textbooks) don't draw much of a distinction between the periodic and aperiodic cases. So I would (probably) still consider Rayleigh's identity to be (more-or-less) Parseval's Theorem, albeit modified to handle periodic functions

"pictures" are pointless

If you want any kind of meaningful response then upload code here using the big green up-arrow in the MaplePrimes toolbar

becuase if you just peel off the sqrt() wrappers in the expression, then you can obtain an answer - see attached.

Hence syntax is not an issue.

With the sqrt() wrappers, *failing* to find a solution would be acceptable (if undesirable) - but it really shouldn't generate an error

restart;
version();
F := x * ( y(x) + x*(x*y(x)) + (x^3*y(x)) );
ode:= diff(y(x),x) = F;
dsolve(ode);

 User Interface: 1298750
         Kernel: 1298750

        Library: 1298750

 

1298750

 

x*(y(x)+x^2*y(x)+x^3*y(x))

 

diff(y(x), x) = x*(y(x)+x^2*y(x)+x^3*y(x))

 

y(x) = _C1*exp((1/20)*x^2*(4*x^3+5*x^2+10))

(1)

 


 

Download odeProb.mw

I have no idea why you want to do what you describe, because it sounds a bit 'circular'. However if(!) I understand the requirement, then maybe it can be met using freeze/thaw, as in


 

restart;
Y:=sqrt(1-y^2)/x;
subExpr:=1-y^2;
Y:=algsubs( subExpr=abs(freeze(subExpr)), Y);
Y:=thaw(%);
Y:=algsubs( abs(subExpr)=subExpr, Y);

(-y^2+1)^(1/2)/x

 

-y^2+1

 

abs(`freeze/R0`)^(1/2)/x

 

abs(-y^2+1)^(1/2)/x

 

(-y^2+1)^(1/2)/x

(1)

 


 

Download freezeThaw.mw

If this rotate+translate process (or some variant of it) produces the actual curves which you wnat, then it would be reasonably simple to incorporate "dummy" x-axis and y-axis along the top and right hand side respectively using some combination of plotttools(line) and textplot()

@michel_h 

You still have pi rather than Pi in the definitions of b[n] and k[n]

SInce the final equation has a single unknown function q(t), but two unknown independent variables 'x' and 't' I still don't think you have any chance of solving this

In you CASE2, the defined function fx is periodic with period [0, T], so in the call to OrthogionalPoly:-FourierSeries(), you need to use x=0..T, rather than x=-T/2..T/2. I have made this simple change in the attached worksheet, and now Parseval's relation seems to hold for your second case.

In both cases, I would suggest that you be a little more careful about the definition of the 'period' of the function.

In CASE1 the function is defined over the range [-1,1] and the Fourier series also uses the range [-1,1] as the basic period. However the associated plot is computes over the range [-2, 2], and you explicitly include two "new" functions defined over the ranges [0,2] and [-2,0]. These "agree" with fx and its Fourier series over the range [-1,1], but not within the ranges [-2,-1] or [1..2]. I find this somewhat confusing. I have 'highlighted' these two definitions in 'red' the attached.

In CASE2, the function fx is defined over the range [0,1] and I have modified the call to FourierSeries() so that it too uses this period. Again you include two "additional" functions in the final plot. One of these is defined over the range [0..1], which is OK, since it is the same as fx over the same range, it is merely superfluous. The second is defined as 0 over the range [-1,0]. Why does this exist? It bears no relation to the periodix function fx, or its Fourier series - It just shouldn't be there! I have highlighted the first of these in 'blue' and the second in 'red' in the attached.

I was't sure what you were trying to do after the plot of your CASE2, so in the attached, I have deleted it


 

restart; with(OrthogonalExpansions)

[BesselSeries, ChebyshevTSeries, ChebyshevUSeries, FourierSeries, GegenbauerSeries, GramSchmidtL2, Haar, HaarSeries, HarmonicWavelet, HarmonicWaveletSeries, HermiteSeries, JacobiSeries, LaguerreSeries, LegendreSeries, MixedSeries, Rational, RationalSeries, RectSeries, SincSeries, SincWavelet, SincWaveletSeries, SphericalSeries, Walsh, WalshSeries, Zernike, ZernikeSeries]

(1)

``

CASE 1:``

fx := piecewise(`and`(-(1/2)*w < x, x < 0), (-1-x/(alpha*tau))*exp(x/(alpha*tau)), 0 < x and x < (1/2)*w, (1-x/(alpha*tau))*exp(-x/(alpha*tau)))

piecewise(-(1/2)*w < x and x < 0, (-1-x/(alpha*tau))*exp(x/(alpha*tau)), 0 < x and x < (1/2)*w, (1-x/(alpha*tau))*exp(-x/(alpha*tau)))

(2)

w := N*alpha*tau; T := M*w

M*N*alpha*tau

(3)

``

a[0] := 0; M := 1; N := 2; alpha := 1; tau := 1; Ck; T; w

2

(4)

Fseries := FourierSeries(fx, x = -(1/2)*w .. (1/2)*w, -n .. n, 'Coefficients'); a := proc (j) options operator, arrow; evalf(eval(op([2, 1, 1], Coefficients), i = j)) end proc; b := proc (j) options operator, arrow; evalf(eval(op([2, 1, 2], Coefficients), i = j)) end proc

proc (j) options operator, arrow; evalf(eval(op([2, 1, 2], Coefficients), i = j)) end proc

(5)

evalf(2*(int(((1-x/(alpha*tau))*exp(-x/(alpha*tau)))^2, x = 0 .. (1/2)*w))); evalf(2*(int(fx^2, x = -(1/2)*w .. (1/2)*w))/w); nterms := [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000, 100000]; seq((1/2)*add(a(k)^2+b(k)^2, k = -p .. p), `in`(p, nterms))

.3938358159, .4125759754, .4243078225, .4282997890, .4303110026, .4315225996, .4319272762, .4321297667, .4322513096, .4322918320, .4323120948, .4323242528, .4323283056

(6)

Plot the results:

tr := [seq(FourierSeries(fx, x = -(1/2)*w .. (1/2)*w, 1 .. k), k = 5 .. 20, 5)]; plots[display](plot([tr[], fx], x = -T .. T, numpoints = 400, axis[1] = [mode = linear], axis[2] = [mode = linear]), plot((1-x/(alpha*tau))*exp(-x/(alpha*tau)), x = 0 .. T, linestyle = dash, thickness = 5, scaling = constrained, color = "Black"), plot(-(1+x/(alpha*tau))*exp(x/(alpha*tau)), x = -T .. 0, linestyle = dash, thickness = 5, scaling = constrained, color = "Black"))

 

Matrix([[k, 1], seq([i, evalf[10](2*Pi*i*(Pi^2*i^2-1+2*(-1)^i*exp(-1))/(Pi^4*i^4+2*Pi^2*i^2+1))], i = 0 .. 1000)])

restart

with(OrthogonalExpansions)

[BesselSeries, ChebyshevTSeries, ChebyshevUSeries, FourierSeries, GegenbauerSeries, GramSchmidtL2, Haar, HaarSeries, HarmonicWavelet, HarmonicWaveletSeries, HermiteSeries, JacobiSeries, LaguerreSeries, LegendreSeries, MixedSeries, Rational, RationalSeries, RectSeries, SincSeries, SincWavelet, SincWaveletSeries, SphericalSeries, Walsh, WalshSeries, Zernike, ZernikeSeries]

(7)

``

CASE 2:``

fx := piecewise(x < 0, 0, 0 < x and x < T, (1-x/(alpha*tau))*exp(-x/(alpha*tau)))

piecewise(x < 0, 0, 0 < x and x < T, (1-x/(alpha*tau))*exp(-x/(alpha*tau)))

(8)

w := N*alpha*tau; T := M*w

M*N*alpha*tau

(9)

``

``

a[0] := 0; M := 1; N := 1; alpha := 1; tau := 1; T; w

1

(10)

Fseries := FourierSeries(fx, x = 0 .. T, -n .. n, 'Coefficients'); a := proc (j) options operator, arrow; evalf(eval(op([2, 1, 1], Coefficients), i = j)) end proc; b := proc (j) options operator, arrow; evalf(eval(op([2, 1, 2], Coefficients), i = j)) end proc

proc (j) options operator, arrow; evalf(eval(op([2, 1, 2], Coefficients), i = j)) end proc

(11)

evalf(2*(int(((1-x/(alpha*tau))*exp(-x/(alpha*tau)))^2, x = 0 .. T))); evalf(2*(int(fx^2, x = 0 .. T))/T); nterms := [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000, 100000]; seq((1/2)*add(a(k)^2+b(k)^2, k = -p .. p), `in`(p, nterms))

.4226913640, .4273910258, .4303260761, .4313241972, .4318270169, .4321299184, .4322310878, .4322817104, .4323120960, .4323222266, .4323272923, .4323303318, .4323313450

(12)

Plot the results:

tr := [seq(FourierSeries(fx, x = 0 .. T, 1 .. k), k = 5 .. 20, 5)]; plots[display](plot([tr[], fx], x = -T .. T, numpoints = 400, axis[1] = [mode = linear], axis[2] = [mode = linear]), plot((1-x/(alpha*tau))*exp(-x/(alpha*tau)), x = 0 .. T, linestyle = dash, thickness = 5, scaling = constrained, color = "Black"), plot(0, x = -T .. 0, linestyle = dash, thickness = 5, scaling = constrained, color = "Black"))

 

NULL


 

Download parseval4.mw

 

Both are simple "rootfinding" conditions, where the "functions" where the desired roots those of s(t)-const, and s(t)-f(s(t)), so that

events=[[s(t)-const], halt]

events=[[s(t)-f(s(t))], halt]

*ought* to work.

A couple of caveats

  1. In both cases there is an implicit assumption, that prior to 'halt' point, s(t) is less than 'const' or f(s(t)) respectively. If you actually want to know (say) the first time that s(t)=const and s(t) is crossing the thrshold 'const' in the positive direction, then you will need to use a conditional trigger, something like events=[[s(t)-const, diff(s(t),t)>0], halt]. In this circumstance, it may or may not be necessary to introduce a dummy variable into you ode system whose value is equal to diff(s(t), t).
  2. It is not obvious (to me) whether you plan on writing f(s(t)) more or less explicitly, or using some function defined defined" externally" to the solution process.

Consider the "toy" example where f(x)=x^2*sin(x). Rather than write

f:=x->x^2*sin(x);
dsolve(odesys, numeric, events=[[s(t)-f(s(t))], halt]

which *might* work (I'm not sure!), I'd probably go with

dsolve(odesys, numeric, events=[[s(t)-s(t)^2*sin(s(t)))], halt]

Why*might* the first version not work? Well, for example, there *might* be an issue with variable "scoping". Consider the whole dsolve() process as running a procedure (written by Maple). As a part of this procedure, the authors *may* have defined a utility function called 'f'. So when the 'events' part of this procedure is run, it will use the local (to the dsolve procedure) definition of 'f', rather than the one you have defined globally. Now you might be able to get around this by calling the 'global' with the :- prefix, as in

f:=x->x^2*sin(x);
dsolve(odesys, numeric, events=[[s(t)-:-f(s(t))], halt]

but it might be safer to avoid this trap (and potentially others like it), by defining f(s(t)) as 'explcitly' as possible - ie the second version given above

If you need further help, it would be better to upload your worksheet (using the big green up-arrow in the Mapleprimes toolbar). It would avoid a lot of 'speculation' on my part!

It would probably be easier if you posted the complete worksheet

The definition of your event puzzles me slightly: it looks like a root-finding event with a conditional trigger.

events = [[[s(t), a*arcsinh(2/a) < s(t)], halt]]

In other words you want to halt when s(t)=0, subject to the condition a*arcsinh(2/a) < s(t)?

So if s(t)=0 ever occurs, you plan to check a*arcsinh(2/a) < s(t), ie a*arcsinh(2/a) <0. I don't think this latter  inequality is ever true, it is positive for all values of 'a', except a=0, when it is 0.

It only makes sense to check this condition if 'a' is a variable - maybe one of the independent variables of the ODE system? Because if 'a' is a parameter, then a*arcsinh(2/a) is a constant, and checking s(t)=0 and s(t)>someConstant would seem to be a bit redundant

you are going to have clarify this (a lot!)

Or provide an example worksheet showing the undesirable behaviour. Worksheets can be uploaded here using the big green up-arrow in the MApleprimes toolbar

  1. I only checked OP's densityplot(), because OP stated "I can not export the .eps file for the densityplot". No mention of problems with other plots
  2. I stated in my original response, that, although I could export the plot as eps - " Resulting plot doesn "look" that great in when loaded into GIMP" So I accept that the result of this process is not a "good" plot
  3. As stated, I used GIMP as the "viewer" and I've always been impressed by the range of different graphics formats which GIMP handles: so if I had to put a bet on it, I'd probably say that the poor quality of the resulting eps plot is caused by the Maple export not the GIMP import - but I can't be 100% certain

I'm running Maple 18 on a 64-bit Windows7 machine. Attached is the eps, of your densityplot which I generated uing the method I described earlier

test.zip

I have had to use a zip file because Mapleprimes does not allow the direct upload of 'eps' files - not sure why:-(

The zip archive just contains one file, test.eps,

 

 

  1. The function fx in your worksheet (defined between -T/2 and T/2) is an odd function, ie f(t)=-f(-t).
  2. cos() is an even function and sin() is an odd function. In a Fourier series representation, all of the cosine term are even, and their sum is even. Similarly all of the sine terms are odd, and their sum is odd.
  3. Since fx is an odd function, its Fourier series representation can only contain odd functions - ie the sine terms. Thus the coefficients of all the cosine terms must be zero. In the conventional notation a[k]*cos()+b[k]*sin(), this means that all the a[k] (except perhaps a[0], which is not multiplied by a 'cosine') must be zero.
  4. If one examines the definition of a[0], this is given by the integral of fx from -T/2 to T/2. Since fx is an odd function, this integral is zero and hence a[0] is zero also

Just noticed that you specified Maple 13. Sorry I can't test something that old, I only keep about the last six versions of Maple on my machine. However it makes me more convinced that it is a version issue!

I can export the densityplot() using the context menu. Right click the plot , then Export->Encapsulated Postscript: this works for me.

Resulting plot doesn "look" that great in when loaded into GIMP - I might consider an alternative format, and I've never been that great at working out which graphics formats are best for which kinds of plots:-(

First 92 93 94 95 96 97 98 Last Page 94 of 207