# Small mismatch in theoretical vs MMA solutions of a damped harmonic oscillator differential equation

I solved the differential equation $$\ddot y + 0.43 \dot y + y =0.1 \cos(2x), \quad y(0)=1,\ \dot y(0)=1 \quad(0)$$ using both the theoretical results from wikipedia and using MMA (I also tested Matlab, but the results are the same of MMA)

DSolve[{y''[x] + 0.43 y'[x] + y[x] == 0.1 Cos[2 x], y'[0] == 1,
y[0] == 1}, y[x], x] // Simplify


The theoretical solution reads $$0.0320427 \cos(3.42077 + 2 x) + 1.59167 \ e^{-0.215 x} \sin(0.70446 + 0.976614 x) \qquad (1)$$ while MMA returns

{{y[x] ->
Cos[0.976614 x] ((1.0308 - 1.77626*10^-19 I) E^(-0.215 x) -
0.0479127 Cos[1.02339 x] + (0.0171106 + 1.77626*10^-19 I) Cos[2.97661 x] +
0.0100658 Sin[1.02339 x] - (0.00123589 + 1.11016*10^-20 I) Sin[2.97661 x]) +
Sin[0.976614 x] ((1.23279 - 1.63692*10^-20 I) E^(-0.215 x) +
0.0100658 Cos[1.02339 x] + (0.00123589 + 1.11016*10^-20 I) Cos[2.97661 x] +
0.0479127 Sin[1.02339 x] + 0.0171106 Sin[2.97661 x])}}


which with a bit of work can be reduced into a simpler formula. In particular, you have to consider the exponential part separately

TrigFactor[
Cos[0.9766140486394818 x] (1.0308020863279805  ) +
Sin[0.9766140486394818 x] (1.2327926138372782 -
1.6369176465928807*^-20 I)  ] E^(-0.215 x)


that gives

(1.60696 - 1.25577*10^-20 I) E^(-0.215 x) Sin[(0.696399 + 0. I) + 0.976614 x]


For the remaining part, with TrigToExp, eliminating the terms lower than 10^18 and with ExpToTrig, in the end it gets

(-0.0320427 + 8.33779*10^-19 I) Sin[(1.84998 + 0. I) + (2. + 0. I) x]


which can be rewritten as

0.0320427 Cos[3.42077 + 2 x]


In the end, MMA solution is

0.0320427 Cos[3.42077 + 2 x] + 1.60696 E^(-0.215 x) Sin[0.696399 + 0.976614 x]


which is slightly different from (1). Is it normal to have this artefact? Does it come from a numerical method to find the solution? I would expect a small difference in the results, but I was expecting a much higher precision in the solution's coefficients.

PS I checked the results also with Matlab, and I get the same solutions of MMA.

EDIT I add the code to find the coefficients of (1) from the theory. Rewrite (0) as $$m \ddot y + c \dot y + k y = Fo \cos( \omega_{dr}x), \quad y(0)=1,\ \dot y(0)=1 \quad(2)$$ that is $$\ddot y + 2\xi \omega \dot y + \omega^2 y = fo \cos( \omega_{dr}x), \quad y(0)=1,\ \dot y(0)=1 \quad(3)$$ with

m = 1;
c = 0.43;
k = 1;
Fo = 0.1;
ωdr = 2;
fo = Fo/m;
ω = Sqrt[k/m];
ξ = c/(2 m ω);
ωd = ω Sqrt[1 - ξ^2];


The solution from the theory is $$y(x)= A e^{-\xi\omega x}\sin(\omega_d x + \theta) + A_0 \cos(\omega_{dr} x - \phi) \quad(4)$$ and have the coefficients evaluated as

    Ao = fo/Sqrt[(ω^2 - ωdr^2)^2 + (2 ξ ω ωdr)^2](*0.0320427*)
ϕ = ArcTan[ω^2 - ωdr^2, 2 ξ ω ωdr] (*2.86241*)
A = Sqrt[(1 + ξ ω - Ao (Sin[ϕ] + ξ ω Cos[ϕ]))^2 + (1 - Ao Cos[ϕ])^2] (*1.59167*)
θ = ArcTan[(1 + ξ ω - Ao (Sin[ϕ] + ξ ω Cos[ϕ])), (1 -  Ao Cos[ϕ])] (*0.70446*)


    Ao = 0.0320427 (*match*)
ϕ = 2.86241 (*match with an appropriate sin/cos conversion*)
A = 1.59167 (*does not match*)
θ = 0.70446 (* does not match *)


The equation to find A and θ has been found imposing the initial condition of eq. (0) on the solution (4). You can find the previous expressions "by hand" or solving the system with MMA

Solve[{A Sin[θ] + Ao Cos[ϕ] == 1, A (-ξ ω) Sin[θ] + A Cos[θ] - Ao Sin[-ϕ] == 1}, {A, θ}]


which gives

{{A -> -1.59167, θ -> -2.43713}, {A -> 1.59167, θ ->
0.70446}}

• I think your Precision (which is machine precision) is not enough, try using arbitrary-precision numbers or exact numbers. – Silvia Feb 5 '15 at 9:48

Your "exact" solution and THE exact solution aren't the same. You probably made some mistakes on your calculations, or a rounding error is lurking in there

realSol = DSolve[{y''[x] + 43/100 y'[x] + y[x] == 1/10 Cos[2 x], y'[0] == 1,  y[0] == 1}, y[x], x]
// Simplify // First


Now, your solution is within 10^-7 of the exact solution, which due to machine precision issues is within 10^-16:

y1[x_] := p.0320427 Cos[3.42077 + 2 x] + 1.60696 E^(-0.215 x) Sin[0.696399 + 0.976614 x]
yy[x_] := Evaluate[y[x] /. realSol]
p1[x_] := y1''[x] + 43/100 y1'[x] + y1[x] - 1/10 Cos[2 x]
pp[x_] := yy''[x] + 43/100 yy'[x] + yy[x] - 1/10 Cos[2 x]

GraphicsRow[{Plot[pp[x], {x, 0, 4 Pi}, PlotRange -> All],
Plot[p1[x], {x, 0, 4 Pi}, PlotRange -> All]}]


• You are saying that realSol is the one the one given by MMA (which match the one from my MMA code), while the solution your[x_] (that I evaluate "by hand") is different from the theoretical one. Plotting the difference do not add anything we already know. – Nicola Feb 6 '15 at 0:38
• So three are the cases, either the theoretical one match the MMA' s or your[x_]  or none of them. If you do the maths, you will see that the coefficient of the exponential (which differs) of your[x_]  match the theoretical one.(see EDIT) – Nicola Feb 6 '15 at 0:45
• @Nicola if you use NDSolve in place of DSolve (for a fully numerical solution), the results will be the same as presented in this answer. The DSolve and NDSolve answers are subtly different (as one would expect), but agree to within $10^{-7}$. Thus, the fact that Mathematica gets the same answer by two different methods, and (if I understood you correctly?) MATLAB gets the same answer again, and substituting Mathematica's symbolic answer back into the equation gives True, strongly suggest that there is a mistake somewhere other than in the software. – Oleksandr R. Feb 6 '15 at 1:55
• I have the same understanding. I posted the "theoretical" method so that someone can checked it and verify it. – Nicola Feb 6 '15 at 2:26