-2
$\begingroup$

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*)

that reads

    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}}
$\endgroup$
  • $\begingroup$ I think your Precision (which is machine precision) is not enough, try using arbitrary-precision numbers or exact numbers. $\endgroup$ – Silvia Feb 5 '15 at 9:48
4
$\begingroup$

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]}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Nicola Feb 6 '15 at 0:38
  • $\begingroup$ 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) $\endgroup$ – Nicola Feb 6 '15 at 0:45
  • $\begingroup$ @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. $\endgroup$ – Oleksandr R. Feb 6 '15 at 1:55
  • $\begingroup$ I have the same understanding. I posted the "theoretical" method so that someone can checked it and verify it. $\endgroup$ – Nicola Feb 6 '15 at 2:26

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.