Part_02.mw

Theory of Oscillations
Part 2: Symbolic Solutions With Initial Conditions

Alexei V. Tikhonenko

Email: avtikhon@okclub.org

General Physics Department, Institute of Nuclear Power Engineering, Russia

This worksheet demonstrates the use of Maple for general symbolic solutions of the equation of oscillation with given initial conditions and how to use ODE solver "dsolve" for these purposes. The general solutions are compared to solutions with initial conditions.

  Introduction.

ODE solver with given initial conditions

To obtain symbolic solutions of the equation of oscillation, we can use the general ODE solver "dsolve" with given initial conditions

Considered situations

1) Free undamped oscillator.

This oscillator is characterized by the frequency w0 and initial conditions.

2) Free damped oscillator.

This oscillator is characterized by the frequency w0, damping factor g and initial conditions.

3) Forced undamped oscillator.

This oscillator is characterized by the frequency w0, amplitude f0 of the external force, the frequency W of the external force and initial conditions.

4) Forced damped oscillator.

This oscillator is characterized by the frequency w0, damping factor g, amplitude f0 of the external force, the frequency W of the external force and the initial conditions.

Initial conditions

Initial conditions for the equation of oscillation are values of the oscillation variable x and its derivative (dx/dt) at the initial moment of time (initial values of the oscillation variable x0 and the velocity v0 of oscillations):

t = 0, x(0) = x0, diff(x(0), t) = v0

For comparison with solutions with initial conditions x_IC(t) and general solution x(t), we will write out x_IC(t) and x(t) for every situation.

  Section 2.1. FREE UNDAMPED OSCILLATOR

Consider the free undamped oscillator. The equation of the oscillation and initial conditions in this situation take the forms:  

diff(x(t), `$`(t, 2))+omega0^2*x(t) = 0

x(0) = x0, diff(x(0), t) = v0

where x(t) is the oscillation variable and w0 is the frequency.

The MAPLE  expression for the differential equation of oscillation, resulting general solution x_01(t) and solution with given initial conditions x_01_IC(t) are:

> restart:
Eq_01:= diff(x(t),t,t)+omega0^2*x(t)=0;
x_01(t):= dsolve(Eq_01);
x_01_IC(t):= dsolve({Eq_01, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_01 := diff(x(t), `$`(t, 2))+omega0^2*x(t) = 0

x_01(t) := x(t) = _C1*sin(omega0*t)+_C2*cos(omega0*t)

x_01_IC(t) := x(t) = v0*sin(omega0*t)/omega0+x0*cos(omega0*t)

The general solution x(t) of the oscillation equation depends on the frequency w0 and two arbitrary constants _C1 and _C2. Solution with given initial conditions x_01_IC(t) depends on the frequency w0, initial value of the oscillation variable x0 and initial velocity v0.

Now we can find the interrelation of the constants _C1 and _C2 and the initial values  of x0 and v0:

> Expr_01:=_C1*sin(omega0*t)+_C2*cos(omega0*t)=1/omega0*v0*sin(omega0*t)+x0*cos(omega0*t);
_C1:=simplify(solve(eval(subs(t=Pi/(2*omega0),Expr_01)),_C1));
_C2:=simplify(solve(eval(subs(t=0,Expr_01)),_C2));

Expr_01 := _C1*sin(omega0*t)+_C2*cos(omega0*t) = v0*sin(omega0*t)/omega0+x0*cos(omega0*t)

_C1 := v0/omega0

_C2 := x0

  Section 2.2. FREE UNDAMPED OSCILLATOR

Consider the free damped oscillator. The equation of oscillation and initial conditions in this situation take the forms:

diff(x(t), `$`(t, 2))+2*gamma*diff(x(t), t)+omega0^2*x(t) = 0

x(0) = x0, diff(x(0), t) = v0

where x(t) is the oscillation variable, w0 is the frequency and g is the damping factor.

The MAPLE  expression for the differential equation of oscillation, resulting general solution x_02(t) and solution with given initial conditions x_02_IC(t) are:

> restart:
Eq_02:= diff(x(t),t,t)+2*gamma*diff(x(t),t)+omega0^2*x(t)=0;
assume(gamma<omega0);
x_02(t):=dsolve(Eq_02);
x_02_IC(t):=dsolve({Eq_02, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_02 := diff(x(t), `$`(t, 2))+2*gamma*diff(x(t), t)+omega0^2*x(t) = 0

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

x_02_IC(t) := x(t) = (x0*gamma+v0)*exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)/(-gamma^2+omega0^2)^(1/2)+x0*exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)

Note that the quantity w0 > g corresponds to the oscillatory solution.

The general solution x(t) of the oscillation equation depends on the frequency w0, damping factor g and two arbitrary constants _C1 and _C2. Solution with given initial conditions x_IC(t) depends on the frequency w0, damping factor g, initial value of oscillation variable x0 and initial velocity v0.

Now we can find the interrelation of the constants _C1 and _C2 and the initial values of x0 and v0:

> Expr_02:=_C1*exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)+_C2*exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)=(x0*gamma+v0)/(-gamma^2+omega0^2)^(1/2)*exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)+x0*exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t);
_C1:=solve(simplify(eval(subs(t=Pi/(sqrt(-gamma^2+omega0^2)),Expr_02)),_C1));
_C2:=simplify(solve(eval(subs(t=0,Expr_02)),_C2));

Expr_02 := ()*exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)+x0*exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t) = (x0*gamma+v0)*exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)/(-gamma^2+omega0^2)^(1/2)+...

_C1 := {omega0 = omega0, x0 = x0}

_C2 := x0

  Section 2.3. FORCED UNDAMPED OSCILLATOR

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

  Sub-Section 2.3.1. FORCED UNDAMPED OSCILLATOR. Nonresonance situation

Consider the forced undamped oscillator with harmonic external force. The equation of oscillation and initial conditions in this situation take the forms:

a)

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

x(0) = x0, diff(x(0), t) = v0

b)

diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_02b*sin(Omega*t)

x(0) = x0, diff(x(0), t) = v0

c)

diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_021c*cos(Omega*t)+f_022c*sin(Omega*t)

x(0) = x0, diff(x(0), t) = v0

d)

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

x(0) = x0, diff(x(0), t) = v0

e)

diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_021e*cos(Omega*t)+f_022e*cos(n*Omega*t)

x(0) = x0, diff(x(0), t) = v0

where x(t) is the oscillation variable, w0 is the frequency, W is the frequency of the external force and a1 and a2 are constant factors.

The MAPLE  expression for the differential equation of oscillation, resulting general solution x(t) and solution with given initial conditions x_IC(t) are:

a)

> restart:
Eq_03a:= diff(x(t),t,t)+omega0^2*x(t)=f_03a*cos(Omega*t);
x_03a(t):=dsolve(Eq_03a);
x_03a_IC(t):=dsolve({Eq_03a, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_03a := diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_03a*cos(Omega*t)

x_03a(t) := x(t) = sin(omega0*t)*_C2+cos(omega0*t)*_C1+f_03a*cos(Omega*t)/(omega0^2-Omega^2)

x_03a_IC(t) := x(t) = sin(omega0*t)*v0/omega0+cos(omega0*t)*(-f_03a+x0*omega0^2-x0*Omega^2)/(omega0^2-Omega^2)+f_03a*cos(Omega*t)/(omega0^2-Omega^2)

Now we can find the interrelation of the constants _C1 and _C2 and the initial values of x0 and v0:

> Expr_03a:=sin(omega0*t)*_C2+cos(omega0*t)*_C1-f_03a*cos(Omega*t)/(-omega0^2+Omega^2)=sin(omega0*t)/omega0*v0+cos(omega0*t)*(x0*Omega^2+f_03a-x0*omega0^2)/(-omega0^2+Omega^2)-f_03a*cos(Omega*t)/(-omega0^2+Omega^2);
_C1:=simplify(solve(eval(subs(t=Pi/omega0,Expr_03a)),_C1));
_C2:=simplify(solve(eval(subs(t=0,Expr_03a)),_C2));

Expr_03a := sin(omega0*t)*_C2+cos(omega0*t)*_C1-f_03a*cos(Omega*t)/(-omega0^2+Omega^2) = sin(omega0*t)*v0/omega0+cos(omega0*t)*(x0*Omega^2+f_03a-x0*omega0^2)/(-omega0^2+Omega^2)-f_03a*cos(Omega*t)/(-o...

_C1 := (-f_03a+x0*omega0^2-x0*Omega^2)/(omega0^2-Omega^2)

_C2 := _C2

b)

> restart:
Eq_03b:= diff(x(t),t,t)+omega0^2*x(t)=f_03b*sin(Omega*t);
x_03b(t):=dsolve(Eq_03b);
x_03b_IC(t):=dsolve({Eq_03b, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_03b := diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_03b*sin(Omega*t)

x_03b(t) := x(t) = sin(omega0*t)*_C2+cos(omega0*t)*_C1+f_03b*sin(Omega*t)/(omega0^2-Omega^2)

x_03b_IC(t) := x(t) = sin(omega0*t)*(-v0*Omega^2-f_03b*Omega+v0*omega0^2)/(omega0*(omega0^2-Omega^2))+cos(omega0*t)*x0+f_03b*sin(Omega*t)/(omega0^2-Omega^2)

Now we can find the interrelation of the constants _C1 and _C2 and the initial values of x0 and v0:

> Expr_03b:=sin(omega0*t)*_C2+cos(omega0*t)*_C1+f_03b*sin(Omega*t)/(omega0^2-Omega^2)=sin(omega0*t)*(-f_03b*Omega+v0*omega0^2-v0*Omega^2)/omega0/(omega0^2-Omega^2)+cos(omega0*t)*x0+f_03b*sin(Omega*t)/(omega0^2-Omega^2);
_C1:=simplify(solve(eval(subs(t=0,Expr_03b)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_03b)),_C2));

Expr_03b := sin(omega0*t)*_C2+cos(omega0*t)*_C1+f_03b*sin(Omega*t)/(omega0^2-Omega^2) = sin(omega0*t)*(-v0*Omega^2-f_03b*Omega+v0*omega0^2)/(omega0*(omega0^2-Omega^2))+cos(omega0*t)*x0+f_03b*sin(Omega...

_C1 := x0

_C2 := _C2

c)

> restart:
Eq_03c:= diff(x(t),t,t)+omega0^2*x(t)=f_031c*cos(Omega*t)+f_032c*sin(Omega*t);
x_03c(t):=dsolve(Eq_03c);
x_03c_IC(t):=dsolve({Eq_03c, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_03c := diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_031c*cos(Omega*t)+f_032c*sin(Omega*t)

x_03c(t) := x(t) = sin(omega0*t)*_C2+cos(omega0*t)*_C1+(-f_031c*cos(Omega*t)-f_032c*sin(Omega*t))/(-omega0^2+Omega^2)

x_03c_IC(t) := x(t) = sin(omega0*t)*(-v0*Omega^2-f_032c*Omega+v0*omega0^2)/(omega0*(omega0^2-Omega^2))+cos(omega0*t)*(-x0*Omega^2-f_031c+x0*omega0^2)/(omega0^2-Omega^2)+(-f_031c*cos(Omega*t)-f_032c*si...

Now we can find the interrelation of the constants _C1 and _C2 and the initial values of x0 and v0:

> Expr_03c:=sin(omega0*t)*_C2+cos(omega0*t)*_C1+(f_031c*cos(Omega*t)+f_032c*sin(Omega*t))/(omega0^2-Omega^2)=sin(omega0*t)*(-v0*Omega^2-f_032c*Omega+v0*omega0^2)/omega0/(omega0^2-Omega^2)+cos(omega0*t)*(-x0*Omega^2-f_031c+x0*omega0^2)/(omega0^2-Omega^2)+(f_031c*cos(Omega*t)+f_032c*sin(Omega*t))/(omega0^2-Omega^2);
_C1:=simplify(solve(eval(subs(t=0,Expr_03c)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_03c)),_C2));

Expr_03c := sin(omega0*t)*_C2+cos(omega0*t)*_C1+(f_031c*cos(Omega*t)+f_032c*sin(Omega*t))/(omega0^2-Omega^2) = sin(omega0*t)*(-v0*Omega^2-f_032c*Omega+v0*omega0^2)/(omega0*(omega0^2-Omega^2))+cos(omeg...

_C1 := (-x0*Omega^2-f_031c+x0*omega0^2)/(omega0^2-Omega^2)

_C2 := _C2

d)

> restart:
Eq_03d:= diff(x(t),t,t)+omega0^2*x(t)=f_03d*cos(n*Omega*t);
x_03d(t):=dsolve(Eq_03d);
x_03d_IC(t):=dsolve({Eq_03d, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_03d := diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_03d*cos(n*Omega*t)

x_03d(t) := x(t) = sin(omega0*t)*_C2+cos(omega0*t)*_C1-f_03d*cos(n*Omega*t)/(-omega0^2+n^2*Omega^2)

x_03d_IC(t) := x(t) = sin(omega0*t)*v0/omega0+cos(omega0*t)*(f_03d-x0*omega0^2+x0*n^2*Omega^2)/(-omega0^2+n^2*Omega^2)-f_03d*cos(n*Omega*t)/(-omega0^2+n^2*Omega^2)

Now we can find the interrelation of the constants _C1 and _C2 and the initial values of x0 and v0:

> Expr_03d:=sin(omega0*t)*_C2+cos(omega0*t)*_C1-f_03d*cos(n*Omega*t)/(-omega0^2+n^2*Omega^2)=sin(omega0*t)/omega0*v0+cos(omega0*t)*(x0*n^2*Omega^2+f_03d-x0*omega0^2)/(-omega0^2+n^2*Omega^2)-f_03d*cos(n*Omega*t)/(-omega0^2+n^2*Omega^2);
_C1:=simplify(solve(eval(subs(t=0,Expr_03d)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_03d)),_C2));

Expr_03d := sin(omega0*t)*_C2+cos(omega0*t)*_C1-f_03d*cos(n*Omega*t)/(-omega0^2+n^2*Omega^2) = sin(omega0*t)*v0/omega0+cos(omega0*t)*(f_03d-x0*omega0^2+x0*n^2*Omega^2)/(-omega0^2+n^2*Omega^2)-f_03d*co...

_C1 := (f_03d-x0*omega0^2+x0*n^2*Omega^2)/(-omega0^2+n^2*Omega^2)

_C2 := _C2

e)

> restart:
Eq_03e:= diff(x(t),t,t)+omega0^2*x(t)=f_031d*cos(Omega*t)+f_032d*cos(n*Omega*t);
x_03e(t):=dsolve(Eq_03e);
x_03e_IC(t):=dsolve({Eq_03e, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_03e := diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_031d*cos(Omega*t)+f_032d*cos(n*Omega*t)

x_03e(t) := x(t) = sin(omega0*t)*_C2+cos(omega0*t)*_C1+((f_032d*Omega^2-f_032d*omega0^2)*cos(n*Omega*t)+(-f_031d*omega0^2+f_031d*n^2*Omega^2)*cos(Omega*t))/(-omega0^4+Omega^2*(n^2+1)*omega0^2-Omega^4*...

x_03e_IC(t) := x(t) = sin(omega0*t)*v0/omega0+cos(omega0*t)*(-f_032d*Omega^2+f_032d*omega0^2+f_031d*omega0^2-f_031d*n^2*Omega^2-x0*omega0^4+x0*omega0^2*n^2*Omega^2+x0*Omega^2*omega0^2-x0*Omega^4*n^2)/...x_03e_IC(t) := x(t) = sin(omega0*t)*v0/omega0+cos(omega0*t)*(-f_032d*Omega^2+f_032d*omega0^2+f_031d*omega0^2-f_031d*n^2*Omega^2-x0*omega0^4+x0*omega0^2*n^2*Omega^2+x0*Omega^2*omega0^2-x0*Omega^4*n^2)/...

Now we can find the interrelation of the constants _C1 and _C2 and the initial values of x0 and v0:

> Expr_03e:=sin(omega0*t)*_C2+cos(omega0*t)*_C1+((f_032d*omega0^2-f_032d*Omega^2)*cos(n*Omega*t)+(f_031d*omega0^2-f_031d*n^2*Omega^2)*cos(Omega*t))/(omega0^4+Omega^2*(-n^2-1)*omega0^2+Omega^4*n^2)=sin(omega0*t)/omega0*v0+cos(omega0*t)*(f_032d*omega0^2-f_032d*Omega^2+f_031d*omega0^2-f_031d*n^2*Omega^2-x0*omega0^4+x0*omega0^2*n^2*Omega^2+x0*Omega^2*omega0^2-x0*Omega^4*n^2)/(-omega0^4+omega0^2*n^2*Omega^2+Omega^2*omega0^2-Omega^4*n^2)+((f_032d*omega0^2-f_032d*Omega^2)*cos(n*Omega*t)+(f_031d*omega0^2-f_031d*n^2*Omega^2)*cos(Omega*t))/(omega0^4+Omega^2*(-n^2-1)*omega0^2+Omega^4*n^2);
_C1:=simplify(solve(eval(subs(t=0,Expr_03e)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_03e)),_C2));

Expr_03e := sin(omega0*t)*_C2+cos(omega0*t)*_C1+((-f_032d*Omega^2+f_032d*omega0^2)*cos(n*Omega*t)+(f_031d*omega0^2-f_031d*n^2*Omega^2)*cos(Omega*t))/(omega0^4+Omega^2*(-n^2-1)*omega0^2+Omega^4*n^2) = ...Expr_03e := sin(omega0*t)*_C2+cos(omega0*t)*_C1+((-f_032d*Omega^2+f_032d*omega0^2)*cos(n*Omega*t)+(f_031d*omega0^2-f_031d*n^2*Omega^2)*cos(Omega*t))/(omega0^4+Omega^2*(-n^2-1)*omega0^2+Omega^4*n^2) = ...

_C1 := (-f_032d*Omega^2+f_032d*omega0^2+f_031d*omega0^2-f_031d*n^2*Omega^2-x0*omega0^4+x0*omega0^2*n^2*Omega^2+x0*Omega^2*omega0^2-x0*Omega^4*n^2)/(-omega0^4+omega0^2*n^2*Omega^2+Omega^2*omega0^2-Omeg...

_C2 := _C2

The general solution x(t) of the oscillation equation depends on the frequency w0, the frequency of external force W and two arbitrary constants _C1 and _C2. Solution with given initial conditions x_IC(t) depends on the frequency w0, the frequency of external force W, initial value of oscillation variable x0 and initial velocity v0.

  Sub-Section 2.3.2. FORCED UNDAMPED OSCILLATOR. Resonance situation

Consider the forced undamped oscillator with harmonic external force of the frequency w0. This is the resonance situation.. The equation of oscillation in this situation takes the forms:

f)

diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_03f*cos(omega0*t)

g)

diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_03g*sin(omega0*t)

h)

diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_031h*cos(omega0*t)+f_032h*sin(omega0*t)

where x(t) is the oscillation variable, w0 is the frequency, W is the frequency of the external force, and a1 and a2 are arbitrary constant factors.

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

f)

> restart:
Eq_03f:= diff(x(t),t,t)+omega0^2*x(t)=f_03f*cos(omega0*t);
x_03f(t):=dsolve(Eq_03f);
x_03f_IC(t):=dsolve({Eq_03f, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_03f := diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_03f*cos(omega0*t)

x_03f(t) := x(t) = sin(omega0*t)*_C2+cos(omega0*t)*_C1+1/2*f_03f*(cos(omega0*t)+sin(omega0*t)*omega0*t)/omega0^2

x_03f_IC(t) := x(t) = sin(omega0*t)*v0/omega0-1/2*cos(omega0*t)*(f_03f-2*x0*omega0^2)/omega0^2+1/2*f_03f*(cos(omega0*t)+sin(omega0*t)*omega0*t)/omega0^2

Now we can find the interrelation of the constants _C1 and _C2 and the initial values of x0 and v0:

> Expr_03f:=sin(omega0*t)*_C2+cos(omega0*t)*_C1+1/2*f_03f*(cos(omega0*t)+sin(omega0*t)*omega0*t)/omega0^2=sin(omega0*t)/omega0*v0+1/2*cos(omega0*t)*(-f_03f+2*x0*omega0^2)/omega0^2+1/2*f_03f*(cos(omega0*t)+sin(omega0*t)*omega0*t)/omega0^2;
_C1:=simplify(solve(eval(subs(t=0,Expr_03f)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_03f)),_C2));

Expr_03f := sin(omega0*t)*_C2+cos(omega0*t)*_C1+1/2*f_03f*(cos(omega0*t)+sin(omega0*t)*omega0*t)/omega0^2 = sin(omega0*t)*v0/omega0+1/2*cos(omega0*t)*(-f_03f+2*x0*omega0^2)/omega0^2+1/2*f_03f*(cos(ome...

_C1 := -1/2*(f_03f-2*x0*omega0^2)/omega0^2

_C2 := _C2

g)

> restart:
Eq_03g:= diff(x(t),t,t)+omega0^2*x(t)=f_03g*sin(omega0*t);
x_03g(t):=dsolve(Eq_03g);
x_03g_IC(t):=dsolve({Eq_03g, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_03g := diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_03g*sin(omega0*t)

x_03g(t) := x(t) = sin(omega0*t)*_C2+cos(omega0*t)*_C1-1/2*f_03g*cos(omega0*t)*t/omega0

x_03g_IC(t) := x(t) = 1/2*sin(omega0*t)*(f_03g+2*v0*omega0)/omega0^2+cos(omega0*t)*x0-1/2*f_03g*cos(omega0*t)*t/omega0

Now we can find the interrelation of the constants _C1 and _C2 and the initial values  of x0 and v0:

> Expr_03g:=sin(omega0*t)*_C2+cos(omega0*t)*_C1-1/2*f_03g*cos(omega0*t)/omega0*t=1/2*sin(omega0*t)*(f_03g+2*v0*omega0)/omega0^2+cos(omega0*t)*x0-1/2*f_03g*cos(omega0*t)/omega0*t;
_C1:=simplify(solve(eval(subs(t=0,Expr_03g)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_03g)),_C2));

Expr_03g := sin(omega0*t)*_C2+cos(omega0*t)*_C1-1/2*f_03g*cos(omega0*t)*t/omega0 = 1/2*sin(omega0*t)*(f_03g+2*v0*omega0)/omega0^2+cos(omega0*t)*x0-1/2*f_03g*cos(omega0*t)*t/omega0

_C1 := x0

_C2 := _C2

h)

> restart:
Eq_03h:= diff(x(t),t,t)+omega0^2*x(t)=f_031h*cos(omega0*t)+f_032h*sin(omega0*t);
x_03h(t):=dsolve(Eq_03h);
x_03h_IC(t):=dsolve({Eq_03h, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_03h := diff(x(t), `$`(t, 2))+omega0^2*x(t) = f_031h*cos(omega0*t)+f_032h*sin(omega0*t)

x_03h(t) := x(t) = sin(omega0*t)*_C2+cos(omega0*t)*_C1+1/2*((f_031h-f_032h*omega0*t)*cos(omega0*t)+sin(omega0*t)*f_031h*omega0*t)/omega0^2

x_03h_IC(t) := x(t) = 1/2*sin(omega0*t)*(f_032h+2*v0*omega0)/omega0^2+1/2*cos(omega0*t)*(-f_031h+2*x0*omega0^2)/omega0^2+1/2*((f_031h-f_032h*omega0*t)*cos(omega0*t)+sin(omega0*t)*f_031h*omega0*t)/omeg...

Now we can find the interrelation of the constants _C1 and _C2 and the initial values  of x0 and v0:

> Expr_03h:=sin(omega0*t)*_C2+cos(omega0*t)*_C1+1/2*((f_031h-f_032h*omega0*t)*cos(omega0*t)+sin(omega0*t)*f_031h*omega0*t)/omega0^2=1/2*sin(omega0*t)*(f_032h+2*v0*omega0)/omega0^2-1/2*cos(omega0*t)/omega0^2*(f_031h-2*x0*omega0^2)+1/2*((f_031h-f_032h*omega0*t)*cos(omega0*t)+sin(omega0*t)*f_031h*omega0*t)/omega0^2;
_C1:=simplify(solve(eval(subs(t=0,Expr_03h)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_03h)),_C2));

Expr_03h := sin(omega0*t)*_C2+cos(omega0*t)*_C1+1/2*((f_031h-f_032h*omega0*t)*cos(omega0*t)+sin(omega0*t)*f_031h*omega0*t)/omega0^2 = 1/2*sin(omega0*t)*(f_032h+2*v0*omega0)/omega0^2-1/2*cos(omega0*t)*...Expr_03h := sin(omega0*t)*_C2+cos(omega0*t)*_C1+1/2*((f_031h-f_032h*omega0*t)*cos(omega0*t)+sin(omega0*t)*f_031h*omega0*t)/omega0^2 = 1/2*sin(omega0*t)*(f_032h+2*v0*omega0)/omega0^2-1/2*cos(omega0*t)*...

_C1 := 1/2*(-f_031h+2*x0*omega0^2)/omega0^2

_C2 := _C2

The general solution x(t) of the oscillation equation depends on the frequency and the frequency of the external force  w0 and two arbitrary constants _C1 and _C2.

  Section 2.4. FORCED DAMPED OSCILLATOR

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

  Sub-Section 2.4.1. FORCED DAMPED OSCILLATOR. Single-frequency external force

Consider the forced damped oscillator with harmonic external force. The equation of oscillation in this situation takes the forms:

a)

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

b)

diff(x(t), `$`(t, 2))+2*gamma*diff(x(t), t)+omega0^2*x(t) = f_04b*sin(sqrt(omega0^2-2*gamma^2)*t)

c)

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

where x(t) is the oscillation variable, w0 is the frequency, W is the frequency of the external force, g is the damping factor and f1, f2 and f3 are constant factors.

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

a)

> restart:
Eq_04a:= diff(x(t),t,t)+2*gamma*diff(x(t),t)+omega0^2*x(t)=f_04a*cos(omega0*t);
assume(gamma<omega0);
x_04a(t):=dsolve(Eq_04a);
x_04a_IC(t):=dsolve({Eq_04a, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_04a := diff(x(t), `$`(t, 2))+2*gamma*diff(x(t), t)+omega0^2*x(t) = f_04a*cos(omega0*t)

x_04a(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+1/2*f_04a*sin(omega0*t)/(gamma*omega0)

x_04a_IC(t) := x(t) = -1/2*exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(f_04a-2*gamma^2*x0-2*v0*gamma)/((-gamma^2+omega0^2)^(1/2)*gamma)+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*x0+1/2*f_04a*...

Now we can find the interrelation of the constants _C1 and _C2 and the initial values  of x0 and v0:

> Expr_04a:=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1+1/2*f_04a/gamma*sin(omega0*t)/omega0=-1/2*exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*(f_04a-2*gamma^2*x0-2*v0*gamma)/(-gamma^2+omega0^2)^(1/2)/gamma+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*x0+1/2*f_04a/gamma*sin(omega0*t)/omega0;
_C1:=simplify(solve(eval(subs(t=0,Expr_04a)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_04a)),_C2));

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

_C1 := x0

_C2 := -1/2*(f_04a-2*gamma^2*x0-2*v0*gamma)/((-gamma^2+omega0^2)^(1/2)*gamma)

b) (Resonance)

> restart:
Eq_04b:= diff(x(t),t,t)+2*gamma*diff(x(t),t)+omega0^2*x(t)=f_04b*cos((sqrt(omega0^2-2*gamma^2))*t);
assume(gamma<omega0);
x_04b(t):=dsolve(Eq_04b);
x_04b_IC(t):=dsolve({Eq_04b, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_04b := diff(x(t), `$`(t, 2))+2*gamma*diff(x(t), t)+omega0^2*x(t) = f_04b*cos((omega0^2-2*gamma^2)^(1/2)*t)

x_04b(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+f_04b*((omega0^2-2*gamma^2)^(1/2)*sin((omega0^2-2*gamma^2)^(1/2)*t)+gamma*cos((...

x_04b_IC(t) := x(t) = -1/2*exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(f_04b-2*x0*gamma^2-2*v0*gamma)/((-gamma^2+omega0^2)^(1/2)*gamma)-1/2*exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*(-2*x0*om...

Now we can find the interrelation of the constants _C1 and _C2 and the initial values  of x0 and v0:

> Expr_04b:=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1-f_04b*(gamma*cos(sqrt(omega0^2-2*gamma^2)*t)+sqrt(omega0^2-2*gamma^2)*sin(sqrt(omega0^2-2*gamma^2)*t))/(-2*gamma*omega0^2+2*gamma^3)=-1/2*exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*(-2*x0*gamma^2-2*v0*gamma+f_04b)/gamma/(-gamma^2+omega0^2)^(1/2)+1/2*exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*(f_04b+2*x0*gamma^2-2*x0*omega0^2)/(gamma^2-omega0^2)-f_04b*(gamma*cos(sqrt(omega0^2-2*gamma^2)*t)+sqrt(omega0^2-2*gamma^2)*sin(sqrt(omega0^2-2*gamma^2)*t))/(-2*gamma*omega0^2+2*gamma^3);
_C1:=simplify(solve(eval(subs(t=0,Expr_04b)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_04b)),_C2));

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

_C1 := -1/2*(-2*x0*omega0^2+f_04b+2*x0*gamma^2)/(-gamma^2+omega0^2)

_C2 := -1/2*(f_04b-2*x0*gamma^2-2*v0*gamma)/((-gamma^2+omega0^2)^(1/2)*gamma)

c)

> restart:
Eq_04c:= diff(x(t),t,t)+2*gamma*diff(x(t),t)+omega0^2*x(t)=f_04c*cos(Omega*t);
assume(gamma<omega0);
x_04c(t):=dsolve(Eq_04c);
x_04c_IC(t):=dsolve({Eq_04c, x(0)=x0, D(x)(0)=v0}, x(t));

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

x_04c(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+f_04c*(2*gamma*Omega*sin(Omega*t)+cos(Omega*t)*omega0^2-Omega^2*cos(Omega*t))/(...

x_04c_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-v0*omega0^4-v0*Omega^4-Omega^4*x0*gamma-4*Omega^2*gamma^3*x0+2*Omega^2*gamma*omega0^2*x0-4*v0*Omega^2*gamma^2-gamma*omega0^4*x0+...x_04c_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-v0*omega0^4-v0*Omega^4-Omega^4*x0*gamma-4*Omega^2*gamma^3*x0+2*Omega^2*gamma*omega0^2*x0-4*v0*Omega^2*gamma^2-gamma*omega0^4*x0+...

Now we can find the interrelation of the constants _C1 and _C2 and the initial values  of x0 and v0:

> Expr_04c:=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1+f_04c*(-Omega^2*cos(Omega*t)+cos(Omega*t)*omega0^2+2*gamma*Omega*sin(Omega*t))/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4)=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*(4*x0*Omega^2*gamma^3-2*gamma*x0*omega0^2*Omega^2+4*v0*Omega^2*gamma^2+gamma*x0*Omega^4-f_04c*Omega^2*gamma+v0*Omega^4+gamma*x0*omega0^4-gamma*f_04c*omega0^2+v0*omega0^4-2*v0*omega0^2*Omega^2)/(-gamma^2+omega0^2)^(1/2)/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4)+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*(-f_04c*omega0^2+f_04c*Omega^2+x0*Omega^4-2*x0*omega0^2*Omega^2+4*x0*Omega^2*gamma^2+x0*omega0^4)/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4)+f_04c*(-Omega^2*cos(Omega*t)+cos(Omega*t)*omega0^2+2*gamma*Omega*sin(Omega*t))/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4);
_C1:=simplify(solve(eval(subs(t=0,Expr_04c)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_04c)),_C2));

Expr_04c := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+f_04c*(2*gamma*Omega*sin(Omega*t)+cos(Omega*t)*omega0^2-Omega^2*cos(Omega*t))/(Omega^4...Expr_04c := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+f_04c*(2*gamma*Omega*sin(Omega*t)+cos(Omega*t)*omega0^2-Omega^2*cos(Omega*t))/(Omega^4...Expr_04c := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+f_04c*(2*gamma*Omega*sin(Omega*t)+cos(Omega*t)*omega0^2-Omega^2*cos(Omega*t))/(Omega^4...

_C1 := (-f_04c*omega0^2+f_04c*Omega^2+x0*Omega^4-2*x0*omega0^2*Omega^2+4*x0*Omega^2*gamma^2+x0*omega0^4)/(omega0^4-2*omega0^2*Omega^2+4*Omega^2*gamma^2+Omega^4)

_C2 := -(-v0*omega0^4-v0*Omega^4-Omega^4*x0*gamma-4*Omega^2*gamma^3*x0+2*Omega^2*gamma*omega0^2*x0-4*v0*Omega^2*gamma^2-gamma*omega0^4*x0+2*v0*omega0^2*Omega^2+f_04c*Omega^2*gamma+gamma*omega0^2*f_04c...

The general solution x(t) of the oscillation equation depends on the frequency w0, damping factor g, and the frequency of the external force, two arbitrary constants _C1 and _C2, the frequency of the external force W and amplitude of the external forces f0.

  Sub-Section 2.4.2. FORCED DAMPED OSCILLATOR. Multifrequency-frequency external force

Consider the forced damped oscillator with harmonic external force. The equation of oscillation in this situation takes the forms:

d)

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

e)

diff(x(t), `$`(t, 2))+2*gamma*diff(x(t), t)+omega0^2*x(t) = f_041e*cos(Omega1*t)+f_042e*cos(Omega2*t)

where x(t) is the oscillation variable, w0 is the frequency, g is damping factor, W, W1, W2 are frequencies of the external forces and f_04d, f_04e are amplitudes of the external forces.

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

d)

> restart:
Eq_04d:=diff(x(t),t,t)+2*gamma*diff(x(t),t)+omega0^2*x(t)=f_041d*cos(Omega*t)+f_042d*cos(2*Omega*t);
assume(gamma<omega0);
x_04d(t):=dsolve(Eq_04d);
x_04d_IC(t):=dsolve({Eq_04d, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_04d := diff(x(t), `$`(t, 2))+2*gamma*diff(x(t), t)+omega0^2*x(t) = f_041d*cos(Omega*t)+f_042d*cos(2*Omega*t)

x_04d(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+1/64*(4*(1/4*Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/4*omega0^4)*(omega0-2*Om...x_04d(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+1/64*(4*(1/4*Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/4*omega0^4)*(omega0-2*Om...x_04d(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+1/64*(4*(1/4*Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/4*omega0^4)*(omega0-2*Om...

x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...x_04d_IC(t) := x(t) = exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*g...

Now we can find the interrelation of the constants _C1 and _C2 and the initial values  of x0 and v0:

> Expr_04d:=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1+1/64*(-16*(Omega+1/2*omega0)*f_042d*(Omega-1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/4*omega0^4)*cos(2*Omega*t)+16*gamma*f_042d*Omega*(1/4*Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/4*omega0^4)*sin(2*Omega*t)+32*f_041d*((-1/2*Omega^2+1/2*omega0^2)*cos(Omega*t)+sin(Omega*t)*Omega*gamma)*(Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/16*omega0^4))/(1/4*Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/4*omega0^4)/(Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/16*omega0^4)=-exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*(-v0*omega0^8-16*v0*Omega^8-20*v0*Omega^2*gamma^2*omega0^4+64*v0*Omega^4*gamma^2*omega0^2-80*v0*Omega^6*gamma^2-64*v0*Omega^4*gamma^4-gamma*omega0^8*x0+gamma*omega0^6*f_042d-16*Omega^8*x0*gamma-80*Omega^6*gamma^3*x0-64*Omega^4*gamma^5*x0+omega0^6*f_041d*gamma+16*f_041d*Omega^6*gamma+16*f_041d*Omega^4*gamma^3+16*gamma^3*Omega^4*f_042d-7*f_041d*Omega^2*gamma*omega0^4+8*f_041d*Omega^4*gamma*omega0^2-7*gamma*Omega^4*f_042d*omega0^2+2*gamma*Omega^2*f_042d*omega0^4+10*v0*omega0^6*Omega^2+40*v0*Omega^6*omega0^2-33*v0*Omega^4*omega0^4+4*gamma*Omega^6*f_042d-33*Omega^4*gamma*omega0^4*x0+64*Omega^4*gamma^3*x0*omega0^2+40*Omega^6*x0*gamma*omega0^2-20*Omega^2*gamma^3*omega0^4*x0+4*Omega^2*gamma^3*f_042d*omega0^2+10*Omega^2*gamma*omega0^6*x0+16*Omega^2*f_041d*gamma^3*omega0^2)/(-gamma^2+omega0^2)^(1/2)/(64*Omega^4*gamma^4+80*Omega^6*gamma^2-64*Omega^4*gamma^2*omega0^2+20*Omega^2*gamma^2*omega0^4+16*Omega^8-40*Omega^6*omega0^2+33*Omega^4*omega0^4-10*omega0^6*Omega^2+omega0^8)+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*(4*f_042d*Omega^6-f_041d*omega0^6-f_042d*omega0^6+16*f_041d*Omega^6+6*f_042d*Omega^2*omega0^4-9*f_042d*omega0^2*Omega^4+9*f_041d*Omega^2*omega0^4-24*f_041d*Omega^4*omega0^2+16*f_041d*Omega^4*gamma^2+16*f_042d*Omega^4*gamma^2-10*x0*omega0^6*Omega^2+33*x0*Omega^4*omega0^4-40*x0*Omega^6*omega0^2+80*x0*Omega^6*gamma^2+64*x0*Omega^4*gamma^4+x0*omega0^8+16*x0*Omega^8-64*x0*Omega^4*gamma^2*omega0^2-16*f_041d*omega0^2*Omega^2*gamma^2-4*f_042d*omega0^2*Omega^2*gamma^2+20*x0*Omega^2*gamma^2*omega0^4)/(64*Omega^4*gamma^4+80*Omega^6*gamma^2-64*Omega^4*gamma^2*omega0^2+20*Omega^2*gamma^2*omega0^4+16*Omega^8-40*Omega^6*omega0^2+33*Omega^4*omega0^4-10*omega0^6*Omega^2+omega0^8)+1/16*(-4*f_042d*(Omega-1/2*omega0)*(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4)*(Omega+1/2*omega0)*cos(2*Omega*t)+4*gamma*Omega*f_042d*(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4)*sin(2*Omega*t)-16*f_041d*((Omega^2-omega0^2)*cos(Omega*t)-2*sin(Omega*t)*Omega*gamma)*(Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/16*omega0^4))/(Omega^4+(-2*omega0^2+4*gamma^2)*Omega^2+omega0^4)/(Omega^4+(-1/2*omega0^2+gamma^2)*Omega^2+1/16*omega0^4);
_C1:=simplify(solve(eval(subs(t=0,Expr_04d)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_04d)),_C2));

Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...Expr_04d := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+1/64*(-16*(Omega-1/2*omega0)*f_042d*(Omega+1/2*omega0)*(1/4*Omega^4+(-1/2*omega0^2+gam...

_C1 := (-4*f_042d*omega0^2*Omega^2*gamma^2-16*f_041d*omega0^2*Omega^2*gamma^2+20*x0*omega0^4*Omega^2*gamma^2-64*x0*omega0^2*Omega^4*gamma^2+16*f_042d*Omega^4*gamma^2+16*f_041d*Omega^4*gamma^2-10*x0*om..._C1 := (-4*f_042d*omega0^2*Omega^2*gamma^2-16*f_041d*omega0^2*Omega^2*gamma^2+20*x0*omega0^4*Omega^2*gamma^2-64*x0*omega0^2*Omega^4*gamma^2+16*f_042d*Omega^4*gamma^2+16*f_041d*Omega^4*gamma^2-10*x0*om..._C1 := (-4*f_042d*omega0^2*Omega^2*gamma^2-16*f_041d*omega0^2*Omega^2*gamma^2+20*x0*omega0^4*Omega^2*gamma^2-64*x0*omega0^2*Omega^4*gamma^2+16*f_042d*Omega^4*gamma^2+16*f_041d*Omega^4*gamma^2-10*x0*om...

_C2 := (-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*gamma*x0*omega0^4*Omega^4-10*gamma*x0*omega0^6*Omega^2-gamma*f_..._C2 := (-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*gamma*x0*omega0^4*Omega^4-10*gamma*x0*omega0^6*Omega^2-gamma*f_..._C2 := (-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*gamma*x0*omega0^4*Omega^4-10*gamma*x0*omega0^6*Omega^2-gamma*f_..._C2 := (-64*gamma^3*x0*omega0^2*Omega^4-4*gamma^3*f_042d*omega0^2*Omega^2+20*gamma^3*x0*omega0^4*Omega^2-40*gamma*x0*omega0^2*Omega^6+33*gamma*x0*omega0^4*Omega^4-10*gamma*x0*omega0^6*Omega^2-gamma*f_...

e)

> restart:
Eq_04e:=diff(x(t),t,t)+2*gamma*diff(x(t),t)+omega0^2*x(t)=f_041e*cos(Omega1*t)+f_042e*cos(Omega2*t);
assume(gamma<omega0);
x_04e(t):=dsolve(Eq_04e);
x_04e_IC(t):=dsolve({Eq_04e, x(0)=x0, D(x)(0)=v0}, x(t));

Eq_04e := diff(x(t), `$`(t, 2))+2*gamma*diff(x(t), t)+omega0^2*x(t) = f_041e*cos(Omega1*t)+f_042e*cos(Omega2*t)

x_04e(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+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(om...x_04e(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+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(om...x_04e(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+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(om...

x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...x_04e_IC(t) := x(t) = -exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omeg...

Now we can find the interrelation of the constants _C1 and _C2 and the initial values of x0 and v0:

> Expr_04e:=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*_C2+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*_C1+((omega0+Omega1)*f_041e*(omega0-Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*cos(Omega1*t)-(Omega1^4+(4*gamma^2-2*omega0^2)*Omega1^2+omega0^4)*f_042e*(Omega2+omega0)*(Omega2-omega0)*cos(Omega2*t)+2*gamma*(Omega1^4+(4*gamma^2-2*omega0^2)*Omega1^2+omega0^4)*f_042e*Omega2*sin(Omega2*t)+2*gamma*Omega1*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*f_041e*sin(Omega1*t))/(Omega1^4+(4*gamma^2-2*omega0^2)*Omega1^2+omega0^4)/(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)=exp(-gamma*t)*sin(sqrt(-gamma^2+omega0^2)*t)*(-gamma*f_042e*omega0^6+4*gamma^3*Omega1^2*x0*Omega2^4+4*gamma^3*x0*Omega2^2*Omega1^4+16*x0*Omega2^2*Omega1^2*gamma^5-4*gamma^3*Omega2^2*f_041e*omega0^2-4*gamma^3*Omega1^2*f_042e*omega0^2+gamma*x0*Omega1^4*omega0^4+4*gamma^3*Omega2^2*x0*omega0^4+2*gamma*f_041e*Omega2^2*omega0^4+4*gamma^3*Omega1^2*x0*omega0^4-2*gamma*x0*omega0^6*Omega2^2+gamma*x0*omega0^4*Omega2^4-gamma*f_041e*Omega2^4*omega0^2+gamma*x0*Omega1^4*Omega2^4-2*gamma*x0*omega0^6*Omega1^2+4*gamma*x0*omega0^4*Omega1^2*Omega2^2+4*v0*Omega2^2*gamma^2*omega0^4+4*v0*omega0^4*Omega1^2*gamma^2-gamma*f_042e*Omega1^4*omega0^2-16*v0*Omega2^2*omega0^2*Omega1^2*gamma^2+16*v0*Omega2^2*gamma^4*Omega1^2+4*v0*Omega2^4*Omega1^2*gamma^2+4*v0*Omega2^2*gamma^2*Omega1^4+2*gamma*f_042e*omega0^4*Omega1^2-gamma*Omega1^2*f_041e*omega0^4-gamma*f_042e*Omega2^2*Omega1^4-gamma*f_042e*Omega2^2*omega0^4-gamma*Omega1^2*f_041e*Omega2^4-4*gamma^3*f_042e*Omega2^2*Omega1^2+v0*omega0^8-4*gamma^3*Omega1^2*f_041e*Omega2^2+2*gamma*f_042e*Omega2^2*omega0^2*Omega1^2+2*gamma*Omega1^2*f_041e*Omega2^2*omega0^2+4*v0*omega0^4*Omega1^2*Omega2^2-2*v0*Omega1^4*Omega2^2*omega0^2-2*v0*omega0^2*Omega1^2*Omega2^4-16*gamma^3*Omega1^2*Omega2^2*omega0^2*x0-2*gamma*x0*Omega1^4*Omega2^2*omega0^2-2*gamma*x0*omega0^2*Omega1^2*Omega2^4-2*v0*omega0^6*Omega1^2+v0*Omega1^4*omega0^4-2*v0*omega0^6*Omega2^2+v0*omega0^4*Omega2^4+v0*Omega1^4*Omega2^4-gamma*f_041e*omega0^6+gamma*x0*omega0^8)/(-gamma^2+omega0^2)^(1/2)/(omega0^4*Omega2^4+4*omega0^4*Omega2^2*gamma^2-2*omega0^6*Omega2^2+omega0^8-2*omega0^2*Omega1^2*Omega2^4-16*omega0^2*Omega1^2*Omega2^2*gamma^2+4*omega0^4*Omega1^2*Omega2^2-2*omega0^6*Omega1^2+4*Omega1^2*gamma^2*Omega2^4+16*Omega1^2*gamma^4*Omega2^2+4*Omega1^2*gamma^2*omega0^4+Omega1^4*Omega2^4+4*Omega1^4*Omega2^2*gamma^2-2*Omega1^4*Omega2^2*omega0^2+Omega1^4*omega0^4)+exp(-gamma*t)*cos(sqrt(-gamma^2+omega0^2)*t)*(x0*omega0^8+f_042e*Omega1^4*Omega2^2-2*f_041e*Omega2^2*omega0^2*Omega1^2-f_042e*Omega1^4*omega0^2+2*f_042e*omega0^4*Omega1^2+f_042e*Omega2^2*omega0^4-f_042e*omega0^6-f_041e*Omega2^4*omega0^2+f_041e*Omega2^4*Omega1^2+2*f_041e*Omega2^2*omega0^4+f_041e*omega0^4*Omega1^2-f_041e*omega0^6-2*f_042e*Omega2^2*omega0^2*Omega1^2-2*x0*omega0^6*Omega1^2-2*x0*Omega1^4*Omega2^2*omega0^2+4*x0*omega0^4*Omega1^2*Omega2^2-2*x0*omega0^2*Omega1^2*Omega2^4+x0*Omega1^4*Omega2^4+x0*omega0^4*Omega2^4-2*x0*omega0^6*Omega2^2+x0*Omega1^4*omega0^4-4*f_041e*Omega2^2*gamma^2*omega0^2+4*f_041e*Omega2^2*gamma^2*Omega1^2+4*f_042e*Omega2^2*gamma^2*Omega1^2-4*f_042e*Omega1^2*gamma^2*omega0^2+4*x0*Omega2^2*gamma^2*Omega1^4+16*x0*Omega2^2*gamma^4*Omega1^2+4*x0*Omega2^2*gamma^2*omega0^4+4*x0*omega0^4*Omega1^2*gamma^2+4*x0*Omega2^4*Omega1^2*gamma^2-16*x0*Omega2^2*omega0^2*Omega1^2*gamma^2)/(omega0^4*Omega2^4+4*omega0^4*Omega2^2*gamma^2-2*omega0^6*Omega2^2+omega0^8-2*omega0^2*Omega1^2*Omega2^4-16*omega0^2*Omega1^2*Omega2^2*gamma^2+4*omega0^4*Omega1^2*Omega2^2-2*omega0^6*Omega1^2+4*Omega1^2*gamma^2*Omega2^4+16*Omega1^2*gamma^4*Omega2^2+4*Omega1^2*gamma^2*omega0^4+Omega1^4*Omega2^4+4*Omega1^4*Omega2^2*gamma^2-2*Omega1^4*Omega2^2*omega0^2+Omega1^4*omega0^4)+1/16*(-4*(1/4*Omega2^4+(-1/2*omega0^2+gamma^2)*Omega2^2+1/4*omega0^4)*f_041e*(omega0+Omega1)*(Omega1-omega0)*cos(Omega1*t)+4*(1/4*Omega1^4+(-1/2*omega0^2+gamma^2)*Omega1^2+1/4*omega0^4)*(Omega2+omega0)*f_042e*(omega0-Omega2)*cos(Omega2*t)+8*gamma*(1/4*Omega1^4+(-1/2*omega0^2+gamma^2)*Omega1^2+1/4*omega0^4)*f_042e*Omega2*sin(Omega2*t)+8*gamma*(1/4*Omega2^4+(-1/2*omega0^2+gamma^2)*Omega2^2+1/4*omega0^4)*sin(Omega1*t)*f_041e*Omega1)/(1/4*Omega2^4+(-1/2*omega0^2+gamma^2)*Omega2^2+1/4*omega0^4)/(1/4*Omega1^4+(-1/2*omega0^2+gamma^2)*Omega1^2+1/4*omega0^4);
_C1:=simplify(solve(eval(subs(t=0,Expr_04e)),_C1));
_C2:=simplify(solve(eval(subs(t=Pi/omega0,Expr_04e)),_C2));

Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...Expr_04e := exp(-gamma*t)*sin((-gamma^2+omega0^2)^(1/2)*t)*_C2+exp(-gamma*t)*cos((-gamma^2+omega0^2)^(1/2)*t)*_C1+(f_041e*(omega0+Omega1)*(Omega2^4+(4*gamma^2-2*omega0^2)*Omega2^2+omega0^4)*(omega0-Om...

_C1 := -(-4*x0*omega0^4*Omega2^2*gamma^2-16*x0*Omega1^2*Omega2^2*gamma^4-4*x0*omega0^4*Omega1^2*gamma^2-4*x0*Omega1^2*Omega2^4*gamma^2-4*x0*Omega1^4*Omega2^2*gamma^2-2*f_041e*Omega2^2*omega0^4-x0*omeg..._C1 := -(-4*x0*omega0^4*Omega2^2*gamma^2-16*x0*Omega1^2*Omega2^2*gamma^4-4*x0*omega0^4*Omega1^2*gamma^2-4*x0*Omega1^2*Omega2^4*gamma^2-4*x0*Omega1^4*Omega2^2*gamma^2-2*f_041e*Omega2^2*omega0^4-x0*omeg..._C1 := -(-4*x0*omega0^4*Omega2^2*gamma^2-16*x0*Omega1^2*Omega2^2*gamma^4-4*x0*omega0^4*Omega1^2*gamma^2-4*x0*Omega1^2*Omega2^4*gamma^2-4*x0*Omega1^4*Omega2^2*gamma^2-2*f_041e*Omega2^2*omega0^4-x0*omeg..._C1 := -(-4*x0*omega0^4*Omega2^2*gamma^2-16*x0*Omega1^2*Omega2^2*gamma^4-4*x0*omega0^4*Omega1^2*gamma^2-4*x0*Omega1^2*Omega2^4*gamma^2-4*x0*Omega1^4*Omega2^2*gamma^2-2*f_041e*Omega2^2*omega0^4-x0*omeg..._C1 := -(-4*x0*omega0^4*Omega2^2*gamma^2-16*x0*Omega1^2*Omega2^2*gamma^4-4*x0*omega0^4*Omega1^2*gamma^2-4*x0*Omega1^2*Omega2^4*gamma^2-4*x0*Omega1^4*Omega2^2*gamma^2-2*f_041e*Omega2^2*omega0^4-x0*omeg...

_C2 := -(-4*gamma^3*Omega1^2*Omega2^4*x0+omega0^6*gamma*f_042e+omega0^6*gamma*f_041e-omega0^8*gamma*x0+2*omega0^6*gamma*Omega2^2*x0-4*omega0^4*gamma^3*Omega1^2*x0-4*omega0^4*gamma^3*Omega2^2*x0-2*omeg...