1
$\begingroup$

I am writing a code that plots an oscillator's response to a periodic excitation. The response consists of the following:

x(t) = xh + xp

where xh is the homogeneous solution and xp is the particular solution. I have already successfully figured out how to plot the particular solution, which is represented as a Fourier series with coefficients a0, an, and bn.

I know that the homogeneous solution is of the form: enter image description here

but I am not sure how to solve for the constants. The initial conditions are that x[0]=0.2 and x'[0]=-0.05 (this is applied for x(t) = xh + xp). I only really know how to solve for these constants when given a differential equation since Mathematica has a function for that, but I am not sure how to ask it to solve for c1 and c2 with both the Fourier series added to the pasted solution above. Everything in the code up to where I define xh correctly works (I apologize for the greek characters in this code I couldn't figure out how to make them look right here, this is my first time posting):

ClearAll["Global`*"]
m = 15; 
k = 960; 
c = 24; 
ωn = Sqrt[k/m]
ζ = c/(2*m*ωn)
ωd = ωn*Sqrt[1 - ζ^2]
Te = 4; 
F = 60; 
F0 = ((F*Te)/4)*t; 
F1 = 90 - 30*t; 
F2 = 0; 
ω = (2*Pi)/Te
r = ω/ωn
MF = 1/((1 - (n*r)^2)^2 + (2*ζ*n*r)^2)^(1/2)
θ = ArcTan[1 - (n*r)^2, 2*ζ*n*r]; 
a0 = (2/Te)*(Integrate[F0, {t, 0, Te/4}] + Integrate[F1, {t, Te/4, (3*Te)/4}] + Integrate[F2, {t, (3*Te)/4, Te}])
an = (2/Te)*(Integrate[F0*Cos[(2*n*Pi*t)/Te], {t, 0, Te/4}] + Integrate[F1*Cos[(2*n*Pi*t)/Te], {t, Te/4, (3*Te)/4}] + Integrate[F2*Cos[(2*n*Pi*t)/Te], {t, (3*Te)/4, Te}])
a[n_] := Simplify[an, Element[n, Integers]]
a[n]
bn = (2/Te)*(Integrate[F0*Sin[(2*n*Pi*t)/Te], {t, 0, Te/4}] + Integrate[F1*Sin[(2*n*Pi*t)/Te], {t, Te/4, (3*Te)/4}] + Integrate[F2*Sin[(2*n*Pi*t)/Te], {t, (3*Te)/4, Te}])
b[n_] := Simplify[bn, Element[n, Integers]]
b[n]
xss[t_] := a0/(2*k) + Sum[(a[n]/k)*MF*Cos[n*ω*t - θ] + (b[n]/k)*Sin[n*ω*t - θ], {n, 1, 100}]
xh[t_] = Exp[(-ζ)*ωn*t]*(c1*Cos[ωd*t] + c2*Sin[ωd*t])
x[t_] = xss[t_] + xh[t_]
Solve[x[0] == 0.2, Derivative[1][x][0] == -0.05, {c1, c2}]
Plot[x[t], {t, 0, 12}]
$\endgroup$
1
$\begingroup$

In general you would use Solve but here it is straightforward:

C1=0.2-xp[0]

and

I let you write the initial equation for C2.

| improve this answer | |
$\endgroup$
1
$\begingroup$
ClearAll["Global`*"]
m = 15;
k = 960;
c = 24;
ωn = Sqrt[k/m];
ζ = c/(2*m*ωn);
ωd = ωn*Sqrt[1 - ζ^2];
Te = 4;
F = 60;
F0 = ((F*Te)/4)*t;
F1 = 90 - 30*t;
F2 = 0;
ω = (2*Pi)/Te;
r = ω/ωn;
MF = 1/((1 - (n*r)^2)^2 + (2*ζ*n*r)^2)^(1/2);
θ = ArcTan[1 - (n*r)^2, 2*ζ*n*r];
a0 = (2/Te)*(Integrate[F0, {t, 0, Te/4}] + 
     Integrate[F1, {t, Te/4, (3*Te)/4}] + 
     Integrate[F2, {t, (3*Te)/4, Te}]);
an = (2/Te)*(Integrate[F0*Cos[(2*n*Pi*t)/Te], {t, 0, Te/4}] + 
     Integrate[F1*Cos[(2*n*Pi*t)/Te], {t, Te/4, (3*Te)/4}] + 
     Integrate[F2*Cos[(2*n*Pi*t)/Te], {t, (3*Te)/4, Te}]);
a[n_] := Simplify[an, Element[n, Integers]];
bn = (2/Te)*(Integrate[F0*Sin[(2*n*Pi*t)/Te], {t, 0, Te/4}] + 
     Integrate[F1*Sin[(2*n*Pi*t)/Te], {t, Te/4, (3*Te)/4}] + 
     Integrate[F2*Sin[(2*n*Pi*t)/Te], {t, (3*Te)/4, Te}]);
b[n_] = Simplify[bn, Element[n, Integers]];
xss[t_] := 
  a0/(2*k) + 
   Sum[(a[n]/k)*MF*Cos[n*ω*t - θ] + (b[n]/k)*
      Sin[n*ω*t - θ], {n, 1, 5}];
xh[t_] := 
  Exp[(-ζ)*ωn*t]*(c1*Cos[ωd*t] + 
     c2*Sin[ωd*t]);
x[t_] := xss[t] + xh[t];
sol = NSolve[{x[0] == 0.2, (D[x[t], t] /. t -> 0) == -0.05}, {c1, c2}];
Plot[x[t] /. sol, {t, 0, 12}]

Mathematica graphics

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.