_

## (dla ambitnych) Przykład w Matlabie

Forum ogólne, ale nie do zadawania "prostych pytań". Wszystko o CAD, czego nie da się bezpośrednio połączyć z tematyką jednego z poniższych forów tematycznych.

### (dla ambitnych) Przykład w Matlabie

Witam wszystkich, których zainteresuje ten temat. Mam do rozwiązania następujące zadanie.
Muszę napisać podobny przykład w Matlabie do tego napisanego w pakiecie Maple.

Kolorem czerwonym zaznaczone są wyniki
Oto on:

Mathieu Functions

The MathieuC and MathieuS Functions
The Mathieu functions MathieuC(a,q,z) and MathieuS(a,q,z) are linearly independent solutions of Mathieu's equation
diff(y(z),(z,2))+(a-2*q*cos(2*z))*y(z);
where a and q are arbitrary parameters. MathieuC and MathieuS are even and odd functions of z. The derivatives of MathieuC and MathieuS with respect to z are MathieuCPrime and MathieuSPrime, respectively. From this definition and the differential equation, the differentiation rules for MathieuCPrime and MathieuSPrime are:
Diff(MathieuCPrime(a,b,f(z)),z)=(diff(f(z),z))*(2*b*cos(2*f(z))-a)*MathieuC(a, b, f(z))

Diff(MathieuSPrime(a,b,f(z)),z)=(diff(f(z),z))*(2*b*cos(2*f(z))-a)*MathieuS(a, b, f(z))
Special values for MathieuC and MathieuS:
[MathieuC(a, b, -z) = MathieuC(a, b, z), MathieuC(a, 0, z) = cos(a^(1/2)*z), MathieuC(a, b, 0) = 1]

[MathieuS(a, b, -z) = -MathieuS(a, b, z), MathieuS(a, 0, z) = sin(a^(1/2)*z)/a^(1/2), MathieuS(a, b, 0) = 0]

For q <> 0 and arbitrary values of the first parameter a, the Mathieu functions are not periodic in z. There exist, however, infinitely many particular values of a, commonly known as characteristic values, such that the Mathieu functions have period k*Pi, where k is a positive integer, corresponding to periodic solutions y(z) = y(k*Pi + z) to Mathieu's equation.
For general complex q <> 0, when a is a characteristic value the Mathieu functions are entire functions of z and the only linearly independent periodic solutions to Mathieu's equation. When the period is Pi; or 2*Pi;, then only one periodic solution exists, otherwise both even and odd periodic solutions exist. In summary, given some (a,q), three different situations can happen:

1. Neither MathieuC(a,q,z) nor MathieuS(a,q,z) is periodic;

2. Mathieu's equation admits a solution with period Pi; or 2*Pi; and only one of MathieuC(a,q,z) or MathieuS(a,q,z) is periodic - the other one is not;

3. Mathieu's equation admits solutions with period 2*n*Pi;, where n >= 2 is a positive integer, and both MathieuC(a,q,z) and MathieuS(a,q,z) are periodic.

The MathieuFloquet Function and the Normalization for MathieuC and MathieuS
According to Floquet's theorem, in all three cases Mathieu's equation admits a solution of the form F(a,q,z) = exp(I*nu(a,q)*z)*P(z);, called a Floquet solution, where nu(a,q); is an arbitrary complex number, commonly known as a characteristic exponent, and P(z); has the same periodicity as the coefficients of Mathieu's equation, namely Pi;. From the invariance of Mathieu's equation under proc (z) options operator, arrow; -z end proc; it follows that F(a,q,-z); is also a solution, which is linearly independent of F(a,q,z); if and only if nu(a,q) is not an integer.
The cases 1, 2, and 3 above respectively correspond to the cases nu(a,q) not a rational number, nu(a,q) an integer, and nu(a,q) a proper rational number, respectively. In case 2, when nu(a,q) is an integer, a is called a characteristic value. Then F(a,q,z); and F(a,q,-z) are proportional and a second (non-periodic) independent solution is constructed as in 20.3.6 and 20.3.7 of reference [1] at the end of this worksheet. For further details, see also Chap. IV in reference [2]. Moreover, the period is Pi; if nu(a,q) is even, and 2*Pi otherwise.
The Maple functions MathieuFloquet(a,q,z) and MathieuFloquetPrime(a,q,z) represent the Floquet solution F(a,q,z) and its derivative with respect to z, respectively. In Maple, MathieuFloquet is normalized so that the L[2]; norm of P(z) is equal to 2*Pi; , and MathieuC and MathieuS are defined in terms of MathieuFloquet, as follows.
> MathieuC(a,q,z) = (MathieuFloquet(a,q,z) + MathieuFloquet(a,q,-z))/(2*MathieuFloquet(a,q,0));
> MathieuS(a,q,z) = (MathieuFloquet(a,q,z) - MathieuFloquet(a,q,-z))/(2*MathieuFloquetPrime(a,q,0));
MathieuC(a, q, z) = 1/2*(MathieuFloquet(a, q, z)+MathieuFloquet(a, q, -z))/MathieuFloquet(a, q, 0)

MathieuS(a, q, z) = 1/2*(MathieuFloquet(a, q, z)-MathieuFloquet(a, q, -z))/MathieuFloquetPrime(a, q, 0)

This yields the following normalization for MathieuC and MathieuS.
> 'MathieuC'(a,q,0) = MathieuC(a,q,0);
> 'MathieuCPrime'(a,q,0) = MathieuCPrime(a,q,0);
> 'MathieuS'(a,q,0) = MathieuS(a,q,0);
> 'MathieuSPrime'(a,q,0) = MathieuSPrime(a,q,0);
MathieuC(a, q, 0) = 1
MathieuCPrime(a, q, 0) = 0
MathieuS(a, q, 0) = 0
MathieuSPrime(a, q, 0) = 1

You have the following special cases when q = 0.
> 'MathieuFloquet'(a,0,z) = MathieuFloquet(a,0,z);
> 'MathieuFloquetPrime'(a,0,z) = MathieuFloquetPrime(a,0,z);

MathieuFloquet(a, 0, z) = exp(a^(1/2)*z*I)
MathieuFloquetPrime(a, 0, z) = a^(1/2)*exp(a^(1/2)*z*I)*I
The Auxiliary Functions MathieuA, MathieuB, and MathieuExponent
Given (a,q), the characteristic exponent nu; entering Floquet solutions is computed using MathieuExponent(a,q). Conversely, given (nu;,q), the values of a entering Floquet solutions and so entering MathieuC(a,q,z) and MathieuS(a,q,z) - regardless of whether nu; is rational and so regardless of whether the functions are periodic - are respectively computed using MathieuA(nu,q) and MathieuB(nu,q). Therefore, the values returned by MathieuA and MathieuB are characteristic values of the corresponding Mathieu function only when nu; is a rational number and MathieuA and MathieuB return a different number only when nu; is an integer (Case 2 in the itemization of section entitled, The MathieuC and MathieuS Functions).
The Mathieu Characteristic functions can be viewed as inverses of each other in that they satisfy: nu = MathieuExponent(MathieuA(nu,q),q) and nu = MathieuExponent(MathieuB(nu,q),q).

The Periodic Solutions MathieuCE and MathieuSE
For a non-negative integer nu;, the even periodic solution corresponding to the nu;th characteristic value is given by MathieuCE(nu,q,z). Similarly, for positive integer nu;, the odd periodic solution corresponding to the nu;th characteristic value is MathieuSE(nu,q,z). Their first derivatives with respect to z are MathieuCEPrime(nu,q,z) and MathieuSEPrime(nu,q,z), respectively. These even and odd periodic Mathieu functions are normalized to have L[2]; norm equal to Pi; if 0 < nu;, and MathieuCE(0,q,z) is normalized to have L[2]; norm equal to 2*Pi;. The following equalities hold.
> MathieuCE(nu,q,z) = MathieuCE(nu,q,0)*MathieuC(MathieuA(nu,q),q,z);
> MathieuSE(nu,q,z) = MathieuSEPrime(nu,q,0)*MathieuS(MathieuB(nu,q),q,z);
MathieuCE(nu, q, z) = MathieuCE(nu, q, 0) MathieuC(MathieuA(nu, q), q, z)
MathieuSE(nu, q, z) = MathieuSEPrime(nu, q, 0) MathieuS(MathieuB(nu, q), q, z)

The even and odd periodic Mathieu functions assume the following special values for q = 0:
> 'MathieuCE'(nu,0,z) = MathieuCE(nu,0,z);
> 'MathieuCEPrime'(nu,0,z) = MathieuCEPrime(nu,0,z);
> 'MathieuSE'(nu,0,z) = MathieuSE(nu,0,z);
> 'MathieuSEPrime'(nu,0,z) = MathieuSEPrime(nu,0,z);
MathieuCE(nu, 0, z) = cos(nu z)
MathieuCEPrime(nu, 0, z) = -nu sin(nu z)
MathieuSE(nu, 0, z) = sin(nu z)
MathieuSEPrime(nu, 0, z) = nu cos(nu z)

Maple recognizes the symmetry and periodicity of MathieuCE and MathieuSE and their derivatives.
> 'MathieuCE'(nu,q,-z) = MathieuCE(nu,q,-z);
> 'MathieuCEPrime'(nu,q,-z) = MathieuCEPrime(nu,q,-z);
> 'MathieuSE'(nu,q,-z) = MathieuSE(nu,q,-z);
> 'MathieuSEPrime'(nu,q,-z) = MathieuSEPrime(nu,q,-z);
MathieuCE(nu, q, -z) = MathieuCE(nu, q, z)
MathieuCEPrime(nu, q, -z) = -MathieuCEPrime(nu, q, z)
MathieuSE(nu, q, -z) = -MathieuSE(nu, q, z)
MathieuSEPrime(nu, q, -z) = MathieuSEPrime(nu, q, z)

> 'MathieuCE'(nu,q,z+2*Pi) = MathieuCE(nu,q,z+2*Pi);
> 'MathieuSE'(nu,q,z+2*Pi) = MathieuSE(nu,q,z+2*Pi);
> MathieuCE(nu,q,z+Pi) assuming nu::even;
> MathieuSE(nu,q,z+Pi) assuming nu::even;
MathieuCE(nu, q, z + 2 Pi) = MathieuCE(nu, q, z)
MathieuSE(nu, q, z + 2 Pi) = MathieuSE(nu, q, z)
MathieuCE(nu, q, z)
MathieuSE(nu, q, z)

A Rational Form for Mathieu's Equation and the Relation with Spheroidal Wave Functions
A rational form of Mathieu's equation, showing the connection between Mathieu and spheroidal wave functions, is obtained by changing variables proc (z) options operator, arrow; arccos(z) end proc;.
> diff(diff(y(z),z),z) = z/(1-z^2)*diff(y(z),z)+((4*z^2-2)*q-a)/(1-z^2)*y(z);

diff(y(z), `\$`(z, 2)) =z*(diff(y(z),z))/(1-z^2)+((4*z^2-2)*q-a)*y(z)/(1-z^2)
Comparing the above with the ODE satisfied by the spheroidal wave functions
> diff(diff(y(z),z),z) = 2*(b+1)*z/(1-z^2)*diff(y(z),z)+((4*q*z^2-c))/(1-z^2)*y(z);
diff(y(z),`\$`(z,2))=2*(b+1)*z*(diff(y(z),z))/(1-z^2)+(4*q*z^2-c)*y(z)/(1-z^2)
you see that Mathieu's equation is the special case of the above when b=-1/2,c=a+2*q. This rational form of Mathieu's equation has two regular singular points at {-1, 1}; and one irregular singular point at infinity;. Therefore, Mathieu functions admit only hypergeometric representation for some particular values of the parameters (a,q).

Series Expansions
Series expansions with respect to the argument are Taylor expansions.
> series(MathieuC(a,q,z), z);

series(1+(q-a/2)*z^2+(-1/3*q+1/6*q^2-1/6*q*a+1/24*a^2)*z^4+O(z^6),z,6)
> series(MathieuS(a,q,z), z=1, 3);

series(MathieuS(a, q, 1)+MathieuSPrime(a, q, 1)*(z-1)+(q*cos(2)-a/2)*MathieuS(a, q, 1)*(z-1)^2+O((z-1)^3),z = 1,3)

> MathieuC(a,q,z) = (MathieuFloquet(a,q,z) + MathieuFloquet(a,q,-z))/(2*MathieuFloquet(a,q,0));
> map(series, %, z, 3);
MathieuC(a, q, z) = 1/2*(MathieuFloquet(a, q, z)+MathieuFloquet(a, q, -z))/MathieuFloquet(a, q, 0)

(series(1+(q-a/2)*z^2+O(z^4),z,4)) = (series(1+(q-a/2)*z^2+O(z^3),z,3))
> series(MathieuCE(a,q,z),z,3);
series(MathieuCE(a, q, 0)+(q-1/2*MathieuA(a, q))*MathieuCE(a, q, 0)*z^2+O(z^3),z,3)
> series(MathieuSEPrime(a,q,z),z,3);
series(MathieuSEPrime(a, q, 0)+(q-1/2*MathieuB(a, q))*MathieuSEPrime(a, q, 0)*z^2+O(z^3),z,3)

Series expansions in the parameter q about 0 and asymptotic expansions are implemented for MathieuCE, MathieuSE, their derivatives, and the characteristic value functions. The series expansions about 0 require the first argument to be an integer, whereas for asymptotic expansions, the first parameter can be symbolic as well.
> series(MathieuCE(1,q,z),q,3);
series(cos(z)-1/8*cos(3*z)*q+(-1/128*cos(z)-1/64*cos(3*z)+1/192*cos(5*z))*q^2+O(q^3),q,3)
> series(MathieuSEPrime(2,q,z),q,3);
series(2*cos(2*z)-1/3*cos(4*z)*q+(-1/144*cos(2*z)+1/64*cos(6*z))*q^2+O(q^3),q,3)
> series(MathieuA(3,q),q);
series(9+1/16*q^2+1/64*q^3+13/20480*q^4-5/16384*q^5+O(q^6),q,6)
> asympt(MathieuSE(1,q,z),q,1);
1/2*sin(z)*exp(-q^(1/2)*cos(z)^2)*HermiteH(0, q^(1/4)*cos(z)*2^(1/2))*Pi^(1/4)*2^(3/4)/(1/q)^(1/8) +Pi^(1/4)*sin(z)*(1/32*HermiteH(0, q^(1/4)*cos(z)*2^(1/2))*2^(3/4)-1/512*2^(3/4)*HermiteH(4, q^(1/4)*cos(z)*2^(1/2))+1/64*HermiteH(2, q^(1/4)*cos(z)*2^(1/2))*2^(3/4))*exp(-q^(1/2)*cos(z)^2)*(1/q)^(3/8) +Pi^(1/4)*sin(z)*(27/2048*HermiteH(0, q^(1/4)*cos(z)*2^(1/2))*2^(3/4)-9/8192*2^(3/4)*HermiteH(4, q^(1/4)*cos(z)*2^(1/2))+5/512*HermiteH(2, q^(1/4)*cos(z)*2^(1/2))*2^(3/4)+1/262144*2^(3/4)*HermiteH(8, q^(1/4)*cos(z)*2^(1/2))-1/16384*HermiteH(6, q^(1/4)*cos(z)*2^(1/2))*2^(3/4))*exp(-q^(1/2)*cos(z)^2)*(1/q)^(7/8) +O(1/q)
> asympt(MathieuB(n,q),q,3);

-2*q+(4*n-2)/(1/q)^(1/2)-(n-1)^2/2-n/2+1/4+
+(3/32*n^2-3/32*n+1/32-1/16*n^3)*(1/q)^(1/2)+(11/256*n-1/16*n^2+
+5/128*n^3-5/256*n^4-3/256)/q+(53/8192-33/4096*n^5-225/8192*n+
+195/4096*n^2-185/4096*n^3+165/8192*n^4)*(1/q)^(3/2)+(189/16384*n^5+
+1413/65536*n-1431/32768*n^2+1575/32768*n^3-2205/65536*n^4-63/16384*n^6-297/65536)/q^2+
+(-6671/262144*n^5-20649/1048576*n+11695/262144*n^2-7613/131072*n^3+48265/1048576*n^4+
+3689/524288*n^6-527/262144*n^7+3961/1048576)*(1/q)^(5/2)+O(1/(q^3))

Plots
The following plots are from reference [1]. Note that the Maple normalization of MathieuCE(0,q,z) differs from [1] by a factor of sqrt(2);.
> plotsetup(inline);
> plot([MathieuA(5-n,q) \$ n=0..5], q=0..15, legend = ['a||(5-n)' \$ n=0..5]);

> plotsetup(inline);
> plot([MathieuB(6-n,q) \$ n=1..5], q=0..15, legend = ['b||(6-n)' \$ n=1..5]);

> plotsetup(inline);
> plot([MathieuCE(n,1,z) \$ n=0..5], z=0..Pi/2, title="Even periodic Mathieu functions, q=1", legend = ['ce||n' \$ n=0..5]);

> plotsetup(inline);
> plot([MathieuSE(n,1,z) \$ n=1..5], z=0..Pi/2, title="Odd periodic Mathieu functions, q=1", legend = ['se||n' \$ n=1..5]);

> plotsetup(inline);
> plot([MathieuCE(n,10,z) \$ n=0..5], z=0..Pi/2, title="Even periodic Mathieu functions, q=10", legend = ['ce||n' \$ n=0..5]);

> plotsetup(inline);
> plot([MathieuSE(n,1,z) \$ n=1..5], z=0..Pi/2, title="Odd periodic Mathieu functions, q=10", legend = ['se||n' \$ n=1..5]);
jurii

Posty: 1
Dołączył(a): lut 13, 2006 18:58