Part_04.mws

Theory of Oscillations
Part 4: Forced Oscillations and Resonance

Alexei V. Tikhonenko

Email: avtikhon@okclub.org

General Physics Department, Institute of Nuclear Power Engineering, Russia

This worksheet demonstrates the use of Maple for analyzing amplitude of forced oscillations.

   Introduction

Amplitudes of forced oscillations

We will analyze amplitudes of forced oscillations for the steady state. In this situation solution x_reg(t) of the equation of oscillation can be presented in forms ( see Section  4.1):

x_reg(t) := A_abs*cos(Omega*t)-A_dis*sin(Omega*t)

or

x_reg(t) := A*cos(Omega*t+phi) = A*cos(Omega*t)*cos(phi)-A*sin(Omega*t)*sin(phi)

where A_abs  and A_dis  are amplitudes of absorption and dispersion, respectively, and A  and j  are the amplitude and phase of the forced oscillations.

Resonance of displacement, velocity and acceleration.

Values A_abs , A_dis , A  and j  are functions of the frequency   of the external force W . We will analyze these functions in detail and find the resonance frequency. In much the same way we will analyze amplitudes of velocity and acceleration.

   Section 4.1. GENERAL SOLUTION  AND ESTABLISHED REGIME

Consider the forced undamped oscillator. Equation of oscillation in this situation takes the form:  

diff(x(t),`$`(t,2))+2*gamma*diff(x(t),t)+omega0^2*x(t) = f0*cos(Omega*t)

where x(t) is the displacement, w 0  is the frequency, g  is   damping factor,   f0  is the amplitude of external force , W  is the frequency   of external force .

The Maple expression for the differential equation of oscillation and the resulting general solution are:

>    restart: interface(showassumed=0):
Eq:= diff(x(t),t,t)+2*gamma*diff(x(t),t)+omega0^2*x(t)=f0*cos(Omega*t);
assume(gamma<omega0);
x_gen(t):=dsolve(Eq);

Eq := diff(x(t),`$`(t,2))+2*gamma*diff(x(t),t)+omega0^2*x(t) = f0*cos(Omega*t)

x_gen(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+((-Omega^2+omega0^2)*cos(Omega*t)+2*gamma*Omega*sin(Omega*t))*f0/(Omega^4+(-2*om...

Note that the quantity w 0  > g  corresponds to the oscillatory solution .

The general solution x_gen(t) is the superposition of the general solution of the homogeneous differential equation x_pr(t) and particular solution of the heterogeneous differential equation x_reg(t):

x_gen(t) = x_pr(t)+x_reg(t)

where

>    x_pr(t):=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1;
x_reg(t):=((-Omega^2+omega0^2)*cos(Omega*t)+2*gamma*Omega*sin(Omega*t))*f0/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4);

x_pr(t) := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1

x_reg(t) := ((-Omega^2+omega0^2)*cos(Omega*t)+2*gamma*Omega*sin(Omega*t))*f0/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4)

x_pr(t) describes the transitory solution, and x_reg(t) describes the steady state:

>    subs([f0=0, Omega=0], x_gen(t));

x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1

that is, x_pr(t) corresponds to the absence of external force;

>    limit(x_pr(t), t=infinity);

0

that is oscillations x_pr(t) damps and system passes to the steady state of oscillations.

   Section 4.2. AMPLITUDES OF ABSORPTION AND AMPLITUDE OF DISPERSION

Consider the steady state of oscillations described by the function:

>    x_reg(t):=((-Omega^2+omega0^2)*cos(Omega*t)+2*gamma*Omega*sin(Omega*t))*f0/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4);

x_reg(t) := ((-Omega^2+omega0^2)*cos(Omega*t)+2*gamma*Omega*sin(Omega*t))*f0/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4)

This function can be represented by the sum:

x_reg(t) := x_S(t)+x_C(t)

where

>    x_S(t):=(2*gamma*Omega*sin(Omega*t))*f0/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4);
x_C(t):=((-Omega^2+omega0^2)*cos(Omega*t))*f0/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4);

x_S(t) := 2*gamma*Omega*sin(Omega*t)*f0/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4)

x_C(t) := (-Omega^2+omega0^2)*cos(Omega*t)*f0/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4)

These functions can be represented by the expressions:

x_reg(t) := A_abs(Omega)*cos(Omega*t)+A_dis(Omega)*sin(Omega*t)

where

>    A_abs(Omega):=2*gamma*Omega*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4);
A_dis(Omega):=(-Omega^2+omega0^2)*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4);

A_abs(Omega) := 2*gamma*Omega*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4)

A_dis(Omega) := (-Omega^2+omega0^2)*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4)

are absorption and dispersion amplitudes of oscillations

   Section 4.3. AMPLITUDE AND PHASE OF  OSCILLATION VARIABLE

Solution x_reg(t) can be represented by the trigonometric expression, too:

x_reg(t) := A*cos(Omega*t+phi)

where A  and j  are the amplitude and phase of forced oscillations.  The last formula can be expanded as follows:

>    x_reg(t):=expand(A(Omega)*cos(Omega*t+phi(Omega)));

x_reg(t) := A(Omega)*cos(Omega*t)*cos(phi(Omega))-A(Omega)*sin(Omega*t)*sin(phi(Omega))

A comparison of formulas

x_reg(t) := A(Omega)*cos(phi(Omega))*cos(Omega*t)-A(Omega)*sin(phi(Omega))*sin(Omega*t)

x_reg(t) := A_abs(Omega)*cos(Omega*t)+A_dis(Omega)*sin(Omega*t)

gives

A_abs(Omega) = A(Omega)*cos(phi(Omega)), A_dis(Omega) = -A(Omega)*sin(phi(Omega))

and

A(Omega)*cos(phi(Omega)) = (-Omega^2+omega0^2)*f0/(Omega^4-2*Omega^2*omega0^2+4*Omega^2*gamma^2+omega0^4), A(Omega)*sin(phi(Omega)) = -2*gamma*Omega*f0/(Omega^4-2*Omega^2*omega0^2+4*Omega^2*gamma^2+ome...

Thus, we have:

>    A(Omega):=simplify(sqrt((A_abs(Omega))^2+(A_dis(Omega))^2));
phi(Omega):=arctan(-A_abs(Omega)/A_dis(Omega));

A(Omega) := (f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))^(1/2)

phi(Omega) := -arctan(2*gamma*Omega/(-Omega^2+omega0^2))

   Section 4.4. AMPLITUDE OF  VELOCITY AND ACCELERATION

To study the amplitudes of velocity and acceleration, we differentiate in t:

>    restart:
x_reg(t):=A(Omega)*cos(Omega*t+phi(Omega));
v_reg(t):=diff(x_reg(t),t);
w_reg(t):=diff(v_reg(t),t);

x_reg(t) := A(Omega)*cos(Omega*t+phi(Omega))

v_reg(t) := -A(Omega)*sin(Omega*t+phi(Omega))*Omega

w_reg(t) := -A(Omega)*cos(Omega*t+phi(Omega))*Omega^2

Therefore amplitudes of velocity and acceleration will have forms:

>    Av(Omega):=A(Omega)*Omega;
Aw(Omega):=A(Omega)*Omega^2;

Av(Omega) := A(Omega)*Omega

Aw(Omega) := A(Omega)*Omega^2

or

>    A(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4));
phi(Omega) := -arctan(2*gamma*Omega/(-Omega^2+omega0^2));
Av(Omega):=A(Omega)*Omega;
Aw(Omega):=A(Omega)*Omega^2;

A(Omega) := (f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))^(1/2)

phi(Omega) := -arctan(2*gamma*Omega/(-Omega^2+omega0^2))

Av(Omega) := (f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))^(1/2)*Omega

Aw(Omega) := (f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))^(1/2)*Omega^2

   Section 4.5. RESONANCE (RESONANCE OF THE DISPLACEMENT )

Now we will analyze variable of oscillations in detail. To do it we will differentiate with respect to W :

>    A(Omega):= sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4));
diff(A(Omega),Omega);

A(Omega) := (f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))^(1/2)

-1/2/(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))^(1/2)*f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4)^2*(-4*omega0^2*Omega+8*Omega*gamma^2+4*Omega^3)

and will solve the result for zero:

>    solve(diff(A(Omega),Omega)=0, Omega);

0, (omega0^2-2*gamma^2)^(1/2), -(omega0^2-2*gamma^2)^(1/2)

Thus function A( W )  has extrema at some points.

We consider two results (since the negative value of the frequency is unsuitable):

>    Omega0:=0;
Omega1:=sqrt(omega0^2-2*gamma^2);

Omega0 := 0

Omega1 := (omega0^2-2*gamma^2)^(1/2)

At these points and at infinity, the function A( W )  is equal to:

>    A(Omega0):=subs(Omega=Omega0,A(Omega));
A(Omega1):=simplify(subs(Omega=Omega1,A(Omega)));
A(infinity):=limit(A(Omega),Omega=infinity);

A(0) := (f0^2/omega0^4)^(1/2)

A((omega0^2-2*gamma^2)^(1/2)) := 1/2/gamma*(f0^2/(omega0^2-gamma^2))^(1/2)

A(infinity) := 0

The function A( W )  is positive and

A(0) < A(sqrt(omega0^2-2*gamma^2))

W1  corresponds to maximum of A( W ) . W1  is called the resonance frequency.

   Section 4.6. RESONANCE OF VELOCITY AND ACCELERATION

Now we will analyze the velocity and acceleration of oscillations in detail. We differentiate with respect to W  :

>    Av(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))*Omega;
Aw(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))*Omega^2;

Av(Omega) := (f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))^(1/2)*Omega

Aw(Omega) := (f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4))^(1/2)*Omega^2

and solve for 0:

>    solve(diff(Av(Omega),Omega)=0, Omega);
solve(diff(Aw(Omega),Omega)=0, Omega);

-omega0, omega0, omega0*I, -I*omega0

0, 1/(omega0^2-2*gamma^2)^(1/2)*omega0^2, -1/(omega0^2-2*gamma^2)^(1/2)*omega0^2

These functions Av( W )  and Av( W )  have extrema:

We consider two results (since negative value of the frequency is unsuitable):

>    Omega2:=omega0;
Omega3:=1/(omega0^2-2*gamma^2)^(1/2)*omega0^2;

Omega2 := omega0

Omega3 := 1/(omega0^2-2*gamma^2)^(1/2)*omega0^2

where W2  is called the resonance frequency of velocity and W3  is called the resonance frequency of acceleration.

   Section 4.7. PLOTS OF AMPLITUDE AND PHASE  FUNCTIONS

The Maple expression for

- frequency w0 ,

- damping factors g ,

- amplitude of external force   f0 ,

are:

>    omega0:=0.75; gamma1:=0.05; gamma2:=0.075; gamma3:=0.1; gamma4:=0.125; gamma5:=0.15; f0:=1;

omega0 := .75

gamma1 := .5e-1

gamma2 := .75e-1

gamma3 := .1

gamma4 := .125

gamma5 := .15

f0 := 1

Amplitude and phase functions for variable of oscillations:

>    A1(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma1^2+Omega^4)):
A2(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma2^2+Omega^4)):
A3(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma3^2+Omega^4)):
A4(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma4^2+Omega^4)):
A5(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma5^2+Omega^4)):
phi1(Omega) := -arctan(2*gamma1*Omega/(-Omega^2+omega0^2)):
phi2(Omega) := -arctan(2*gamma2*Omega/(-Omega^2+omega0^2)):
phi3(Omega) := -arctan(2*gamma3*Omega/(-Omega^2+omega0^2)):
phi4(Omega) := -arctan(2*gamma4*Omega/(-Omega^2+omega0^2)):
phi5(Omega) := -arctan(2*gamma5*Omega/(-Omega^2+omega0^2)):
A_abs1(Omega):=2*gamma1*Omega*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma1^2+Omega^4):
A_abs2(Omega):=2*gamma2*Omega*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma2^2+Omega^4):
A_abs3(Omega):=2*gamma3*Omega*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma3^2+Omega^4):
A_abs4(Omega):=2*gamma4*Omega*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma4^2+Omega^4):
A_abs5(Omega):=2*gamma5*Omega*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma5^2+Omega^4):

A_dis1(Omega):=(-Omega^2+omega0^2)*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma1^2+Omega^4):
A_dis2(Omega):=(-Omega^2+omega0^2)*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma2^2+Omega^4):
A_dis3(Omega):=(-Omega^2+omega0^2)*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma3^2+Omega^4):
A_dis4(Omega):=(-Omega^2+omega0^2)*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma4^2+Omega^4):
A_dis5(Omega):=(-Omega^2+omega0^2)*f0/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma5^2+Omega^4):

plot([A1(Omega),A2(Omega),A3(Omega),A4(Omega),A5(Omega)],Omega=0..1.6,color=[red,blue,magenta,gold, navy],linestyle=[1,3,3,3,1],thickness=[2,1,1,1,2],axes=FRAME);
plot([phi1(Omega),phi2(Omega),phi3(Omega),phi4(Omega),phi5(Omega)],Omega=0..1.6,color=[red,blue,magenta,gold,navy],linestyle=[1,3,3,3,1],thickness=[2,1,1,1,2],axes=FRAME);
plot([A_abs1(Omega),A_abs2(Omega),A_abs3(Omega),A_abs4(Omega),A_abs5(Omega)],Omega=0..1.6,color=[red,blue,magenta,gold,navy],linestyle=[1,3,3,3,1],thickness=[2,1,1,1,2],axes=FRAME);
plot([A_dis1(Omega),A_dis2(Omega),A_dis3(Omega),A_dis4(Omega),A_dis5(Omega)],Omega=0..1.6,color=[red,blue,magenta,gold,navy],linestyle=[1,3,3,3,1],thickness=[2,1,1,1,2],axes=FRAME);

[Maple Plot]

[Maple Plot]

[Maple Plot]

[Maple Plot]

   Section 4.8. PLOTS OF AMPLITUDE FUNCTIONS FOR VELOCITY AND ACCELERATION

The Maple expression for

- frequency w0 ,

- damping factors g ,

- amplitude of external force   f0 ,

are:

>    omega0:=0.75; gamma1:=0.05; gamma2:=0.075; gamma3:=0.1; gamma4:=0.125; gamma5:=0.15; f0:=1;

omega0 := .75

gamma1 := .5e-1

gamma2 := .75e-1

gamma3 := .1

gamma4 := .125

gamma5 := .15

f0 := 1

Amplitude functions for variable of oscillations, velocity and acceleration:

>    Av1(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma1^2+Omega^4))*Omega:
Av2(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma2^2+Omega^4))*Omega:
Av3(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma3^2+Omega^4))*Omega:
Av4(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma4^2+Omega^4))*Omega:
Av5(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma5^2+Omega^4))*Omega:
Aw1(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma1^2+Omega^4))*Omega^2:
Aw2(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma2^2+Omega^4))*Omega^2:
Aw3(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma3^2+Omega^4))*Omega^2:
Aw4(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma4^2+Omega^4))*Omega^2:
Aw5(Omega):=sqrt(f0^2/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma5^2+Omega^4))*Omega^2:
plot([A1(Omega),A2(Omega),A3(Omega),A4(Omega),A5(Omega)],Omega=0..1.6,color=[red,blue,magenta,gold, navy],linestyle=[1,3,3,3,1],thickness=[2,1,1,1,2],axes=FRAME);
plot([Av1(Omega),Av2(Omega),Av3(Omega),Av4(Omega),Av5(Omega)],Omega=0..1.6,color=[red,blue,magenta,gold, navy],linestyle=[1,3,3,3,1],thickness=[2,1,1,1,2],axes=FRAME);
plot([Aw1(Omega),Aw2(Omega),Aw3(Omega),Aw4(Omega),Aw5(Omega)],Omega=0..1.6,color=[red,blue,magenta,gold, navy],linestyle=[1,3,3,3,1],thickness=[2,1,1,1,2],axes=FRAME);

[Maple Plot]

[Maple Plot]

[Maple Plot]

>   

In this Section we consider two situations: Resonance oscillations and Nonresonance oscillations.

   Section 4.9. PLOTS OF DISPLACEMENT, VELOCITY AND ACCELERATION

The Maple expressionS for

- frequency w0 ,

- damping factors g ,

- amplitude of external force   f0 ,
- frequencies of external force
W,

- initial conditions and constants

are:

>    restart:
omega0:=1.75; gamma:=0.15; f0:=1; Omega01:=1.25; Omega02:=sqrt(omega0^2-2*gamma^2); Omega03:=2;
x0:=0; v0:=0;
_C1 := x0;
_C2 := -1/2*(f0-2*gamma^2*x0-2*v0*gamma)/(-gamma^2+omega0^2)^(1/2)/gamma;

omega0 := 1.75

gamma := .15

f0 := 1

Omega01 := 1.25

Omega02 := 1.737095277

Omega03 := 2

x0 := 0

v0 := 0

_C1 := 0

_C2 := -1.911797782

Plots of displacement, velocity and acceleration for different frequencies of external force ( W02  is the resonance frequency):

>    x_gen1(t):=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1+f0*(cos(Omega01*t)*omega0^2+2*gamma*Omega01*sin(Omega01*t)-Omega01^2*cos(Omega01*t))/(Omega01^4+(-2*omega0^2+4*gamma^2)*Omega01^2+omega0^4):
x_gen2(t):=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1+f0*(cos(Omega02*t)*omega0^2+2*gamma*Omega02*sin(Omega02*t)-Omega02^2*cos(Omega02*t))/(Omega02^4+(-2*omega0^2+4*gamma^2)*Omega02^2+omega0^4):
x_gen3(t):=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1+f0*(cos(Omega03*t)*omega0^2+2*gamma*Omega03*sin(Omega03*t)-Omega03^2*cos(Omega03*t))/(Omega03^4+(-2*omega0^2+4*gamma^2)*Omega03^2+omega0^4):

v_gen1(t):=diff(x_gen1(t),t):
v_gen2(t):=diff(x_gen2(t),t):
v_gen3(t):=diff(x_gen3(t),t):

w_gen1(t):=diff(v_gen1(t),t):
w_gen2(t):=diff(v_gen2(t),t):
w_gen3(t):=diff(v_gen3(t),t):

plot([x_gen1(t),x_gen2(t),x_gen3(t)],t=0..60,color=[red,blue,magenta],linestyle=[1,3,3],thickness=[2,1,1],axes=FRAME);
plot([v_gen1(t),v_gen2(t),v_gen3(t)],t=0..60,color=[red,blue,magenta],linestyle=[1,3,3],thickness=[2,1,1],axes=FRAME);
plot([w_gen1(t),w_gen2(t),w_gen3(t)],t=0..60,color=[red,blue,magenta],linestyle=[1,3,3],thickness=[2,1,1],axes=FRAME);

[Maple Plot]

[Maple Plot]

[Maple Plot]

© 2002. Alexei V. Tikhonenko