Message Boards Message Boards

Solve problem with a particular initial condition using NDSolve?

Posted 8 years ago

.enter image description here

I attached my program. And the picture above is the system .It's a spinning disk on a movable car, and give it a base-excitation. I'm using the 1-2-3 Euler angle to describe the configuration of the disk, so alpha=beta=0 means the rod between disk and car is completely straight.

I want to figure out that whether the increase of the spin rate can reduce the vibration of the car or not. But when I set initial condition alpha=beta=0 the program pop up an error message.

I've check the determinant [M] , it's not zero, so it shouldn't be a singular point.

Maybe setting alpha[0] and beta[0] equal to 0.0001 or a even smaller number is a easy way, but I still want to find out the answer of this situation.

Attachments:
POSTED BY: Homer Chiou
4 Replies
Posted 8 years ago

Thanks for your reply again, I found out that the term 2 L*Yc'[t] ([Beta]'[t] Sin[[Alpha][t]] Sin[[Beta][t]] - [Alpha]'[ t] Cos[[Alpha][t]] Cos[[Beta][t]])) causes the problem. So when mathematica reduce this into normal form, turns it into Cot[[Alpha][0]] and Cot[[Beta][0]]? Thanks again ^^

POSTED BY: Homer Chiou

You can copy the code, paste it into a notebook and evaluate. The first, long part is simply the differential system (copied and pasted from your notebook), where I have just separated the equations from the initial conditions. The first thing NDSolve does it to reduce the equations into normal form, by solving for the second derivatives. This is where the problem arises. The expression that Mathematica gets contains terms, that are singular precisely in your initial data. What I suggested is to solve for the second derivatives in another way, that does not produce the problematic terms. Then we can feed the safe expression into NDSolve, which does its job.

POSTED BY: Gianluca Gorni

Your determinant is nonzero but rather small (10^-7). However, it seems that the problem is symbolic. When Mathematica solves for the second derivatives, unfortunately the expression contains Cot[\[Alpha][0]] and Cot[\[Beta][0]], which are infinite at your initial conditions. One way to avoid the problem is to manually invert the matrix and thus bypass the built-in solver:

eqDiff = {0.098` Cos[\[Beta][t]] Sin[\[Alpha][t]] - \[Alpha][t] - 
     0.0014074074074074073` Sin[2 \[Beta][t]] Derivative[1][\[Alpha]][
       t] Derivative[1][\[Beta]][t] - 
     1/200 Cos[\[Beta][t]] Derivative[1][\[Beta]][t] Derivative[
       1][\[Phi]][t] + 
     0.010000000000000002` Cos[\[Alpha][t]] Cos[\[Beta][t]] (
       Yc^\[Prime]\[Prime])[t] - 
     0.0035925925925925925` Cos[\[Beta][
        t]]^2 (\[Alpha]^\[Prime]\[Prime])[t] - 
     1/200 Sin[\[Beta][t]]^2 (\[Alpha]^\[Prime]\[Prime])[t] - 
     1/200 Sin[\[Beta][t]] (\[Phi]^\[Prime]\[Prime])[t] == 0, 
   0.098` Cos[\[Alpha][t]] Sin[\[Beta][t]] - \[Beta][t] + 
     0.0014074074074074073` Cos[\[Beta][t]] Sin[\[Beta][
        t]] Derivative[1][\[Alpha]][t]^2 + 
     1/200 Cos[\[Beta][t]] Derivative[1][\[Alpha]][t] Derivative[
       1][\[Phi]][t] - 
     0.010000000000000002` Cos[\[Beta][t]] (Xc^\[Prime]\[Prime])[t] - 
     0.010000000000000002` Sin[\[Alpha][t]] Sin[\[Beta][t]] (
       Yc^\[Prime]\[Prime])[t] - 
     0.0035925925925925925` (\[Beta]^\[Prime]\[Prime])[t] == 0, 
   1/200 (-Cos[\[Beta][t]] Derivative[1][\[Alpha]][t] Derivative[
         1][\[Beta]][t] - 
       Sin[\[Beta][t]] (\[Alpha]^\[Prime]\[Prime])[
         t] - (\[Phi]^\[Prime]\[Prime])[t]) == 0, 
   3.` + 10 (t UnitStep[t] - 
        Sin[t] UnitStep[t] + (2 \[Pi] - t) UnitStep[-2 \[Pi] + t] + 
        Sin[t] UnitStep[-2 \[Pi] + t]) - 10 Xc[t] + 
     0.010000000000000002` Sin[\[Beta][t]] Derivative[1][\[Beta]][
       t]^2 - 1.1` (Xc^\[Prime]\[Prime])[t] - 
     0.010000000000000002` Cos[\[Beta][t]] (\[Beta]^\[Prime]\[Prime])[
       t] == 0, -Yc[t] - (Yc^\[Prime]\[Prime])[t] - 
     0.1` (0.1` Cos[\[Beta][t]] Sin[\[Alpha][t]] Derivative[
          1][\[Alpha]][t]^2 + 
        0.2` Cos[\[Alpha][t]] Sin[\[Beta][t]] Derivative[1][\[Alpha]][
          t] Derivative[1][\[Beta]][t] + 
        0.1` Cos[\[Beta][t]] Sin[\[Alpha][t]] Derivative[1][\[Beta]][
          t]^2 + (Yc^\[Prime]\[Prime])[t] - 
        0.1` Cos[\[Alpha][t]] Cos[\[Beta][
           t]] (\[Alpha]^\[Prime]\[Prime])[t] + 
        0.1` Sin[\[Alpha][t]] Sin[\[Beta][
           t]] (\[Beta]^\[Prime]\[Prime])[t]) == 0};
inCond = {Xc[0] == 0.3`, Derivative[1][Xc][0] == 0, \[Beta][0] == 0, 
   Derivative[1][\[Beta]][0] == 0, \[Phi][0] == 0, 
   Derivative[1][\[Phi]][0] == 200, \[Alpha][0] == 0, 
   Derivative[1][\[Alpha]][0] == 0, Yc[0] == 0, 
   Derivative[1][Yc][0] == 0};
vars = {\[Alpha], \[Beta], \[Phi], Yc, Xc};
vars2 = D[Through[vars[t]], t, t];
m = Map[Coefficient[Map[First, eqDiff], #] &, vars2];
b = Map[First, eqDiff] /. Thread[vars2 -> 0];
linSol1 = Inverse[m].(-b);
NDSolve[{vars2 == linSol1, inCond}, vars, {t, 0, 100}]

The built-in linear solvers have the problem:

linSol2 = LinearSolve[m, -b];
linSol3 = Solve[eqDiff, vars2];
{linSol1, linSol2, linSol3} /. t -> 0 /. Solve[inCond]
POSTED BY: Gianluca Gorni
Posted 8 years ago

Thanks for your reply, but the code is a little bit hard to read. Do you mind upload the mathematica file or send it to me ? Thanks alot. I'm a little confused that the initial condition makes Cot[[Alpha][0]] and Cot[[Beta][0]] infinite results this problem.Please explain a little more for me to understand, I need to understand more detail of this problem, thanks again.

POSTED BY: Homer Chiou
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract