dsolve isn't solving my symbolic ODE correctly

I'm trying so solve two simple nonhomogeneous ODEs using dsolve. My forcing functions always have the some form, but I am allowing some parameters of the functions to change, i.e. the frequencies and the amplitude. The forced ODEs are:

and

I expect the particular solutions to be sums of sinusoids at the forcing frequencies. If we just take the Fourier transform of the first equation, we find that:

which is just an amplitude gain with no phase shift. The symbolic toolbox sometimes returns a good particular solution (coefficients are actual numbers in my code):

beta_s_sol(t) =
 
- 0.0444*cos(2.0944*t) - 0.0934*cos(3.1416*t)

but other times it's not correct. For the same system parameters, the gamma solution was:

gamma_s_sol(t) =
 
-7.7215e-11*exp(-8.6936*t)*(5.7469e+08*exp(8.6936*t)*sin(2.0944*t) - 1.2101e+09*exp(8.6936*t)*sin(3.1416*t))

Is there a known bug in the symbolic toolbox? I'm designing a control system that relies on those particular solutions and when I get two good solutions, the control system works very well. When I don't get the right solution, the control system doesn't work. I may have some fundamental issues elsewhere, but the controller response correlating with the ODE solution issue is pretty convincing to me.

[Torsten]

Is a/b resp. c/d above always positive ? Otherwise, you'll get exp-terms in the solution of the homogenous system.

Answer 1

[John D'Errico]

Look carefully at your solution, because you missed something.

gamma_s_sol(t) =
-7.7215e-11*exp(-8.6936*t)*(5.7469e+08*exp(8.6936*t)*sin(2.0944*t) - 1.2101e+09*exp(8.6936*t)*sin(3.1416*t))
            ______________             _____________

So we see a product of terms. In one case, there is a negative exponential. In the second, we have a positive exponential.

Just use simplify. I also used vpa, to turn the fractions into numbers, making it easier to read. And your numbers were only given to a few decimals anyway, so 5 digits is about all those expressions are worth.

vpa(simplify(gamma_s_sol),5)
ans(t) =
0.093438*sin(3.1416*t) - 0.044375*sin(2.0944*t)

Just a sum of simple trig terms. The exponentials cancel out. All you needed to do was use simplify. The exponentials looked like they were there. But really, they were not, a pure illusion done with mirrors. It was easy to miss. And no bug. Along the way, I assume that somewhere along the way, something exponential gets generated, that ends up completely cancelling out for your problem.

[Dominic Riccoboni]

Thank you for looking into this in detail. I agree that doing so for the above gamma result gets a good solution. I did actually try simplify at one point and the exponential terms went away, I must have just forgotten to do that before writing this post. It probably would have prevented the post in the first place.

I think I was just frustrated – there’s got to be some other factor I’m changing too readily that’s causing the control system to break and I keep trying to blame the wonky looking solutions I sometimes get when I don’t use simplify. I’ll just have simplify in there permanently and that’ll probably be that.

Answer 2

[Torsten]

Is a/b resp. c/d above always positive ? Otherwise, you'll get exp-terms in the solution of the homogenous system.

Maybe you can just copy the solution here and integrate it as-is in your program code:

https://www.wolframalpha.com/input?i=a*y%27%27%28x%29+%2B+b*y%28x%29+%3D+d*%28-w%5E2*Sin%5Bw*x%5D%2Bv%5E2*Sin%5Bv*x%5D%29

[Dominic Riccoboni]

Thank you for your response Torsten. I'm actually setting the homogeneous response to zero by zeroing out the arbitrary constants C1 and C2 in the result.

[Torsten]

And why do you make permanent calls to "dsolve" if you could use the above analytical solution for your 2nd order ODE as fixed part of your code ?

[Dominic Riccoboni]

The forcing function I showed here was just a concrete example. My code currently solves that ODE for arbitrary forcing functions (well, within reason).

[Torsten]

And these forcing functions do not allow a symbolic solution without giving numerical values to the parameters involved ?

[Dominic Riccoboni]

No, they definitely could, it's just the way I've structured my code, I don't want someone to have to select from a list of parameterized solutions. I just want them to be able to ask for a desired trajectory in x and y coordinates (the forcing functions are the x and y accelerations) and the code will automatically compute consistent trajectories in the angles (gamma and beta). It's just to be as general as possible.