Circular reference in differential equations

Questions
, , ,
1

Hi Guys

Please refer to the below code. Each differential equation presents acceleration terms from the other 2 differential equations, I assume that this would be a circular reference. I cannot get the system to solve, can anyone please guide me as to what is wrong?

from wolframclient.evaluation import WolframLanguageSession
from wolframclient.language import wl, wlexpr
import matplotlib.pyplot as plt
import numpy as np

session = WolframLanguageSession()

wolfram_code = """

                    IC = {x[0] == 0, theta[0] == 0, phi[0] == 0, x'[0] == 0, theta'[0] == 0, phi'[0] == 0};

                    


                    
                    l = 1;
                    m = 1000;
                    r = 0.1;
                    g = 9.81;

                    Kt = 1000;
                    Kb = 10000;
                    Ib = 1000;
                    It = 1000;

                    C1 = 0;
                    C2 = 0;
                    C3 = 0;



                    one = -0.5*m*(2*l*theta''[t]*Cos[theta[t]]-2*l*theta'[t]*Sin[theta[t]]+2*r*phi''[t]*Sin[phi[t]]+2*r*phi'[t]^2*Cos[theta[t]])-C1*x'[t];


                    two_part_one = -0.5*m*(2*l*x''[t]*Cos[theta[t]]-2*l*x'[t]*theta'[t]*Sin[theta[t]]+2*r*phi''[t]*Cos[theta[t]]*Sin[theta[t]]+2*r*l*theta'[t]*phi'[t]*Cos[2*theta[t]]);
                    two_part_two = 0.5*m*(-2*l*theta'[t]*x'[t]*Sin[theta[t]]+2*r*l*theta'[t]*phi'[t]*Cos[2*theta[t]])+Kb*theta[t]+m*g*l*Sin[theta[t]]-C2*theta'[t];
                    two = two_part_one+two_part_two;

                    three_part_one = -0.5*m*(m*r*x''[t]*Sin[theta[t]]+2*r*x'[t]*theta'[t]*Cos[theta[t]]+2*r*l*theta''[t]*Cos[theta[t]]*Sin[theta[t]]+2*r*l*theta'[t]^2*Cos[theta[t]]);
                    three_part_two = 0.5*m*(-2*r*phi'[t]*x'[t]*Cos[phi[t]])-Kt*phi[t]-C3*phi'[t];
                    three = three_part_one+three_part_two; 

                    S = NDSolve[{x''[t] == one/m, theta''[t] == two/(m*l^2+2*Ib), phi''[t] == three/(m*r+2*It),IC},{x,theta,phi},{t,0,10}]



"""

# Evaluate the Wolfram code
result = session.evaluate(wlexpr(wolfram_code))

# Terminate the Wolfram session
session.terminate()
2

The below is the error which I received.

Tag Times in Times[Blank[Times[-500.0, Plus[Times[0.2, Cos[theta[t]], Power[Derivative[1][phi][t], 
2]], Times[-2, Sin[theta[t]], Derivative[1][theta][t]], Times[0.2, Sin[phi[t]], Derivative[2][phi][t]], Times[2, Cos[theta[t]], Derivative[2][theta][t]]]]], Pattern[two, Blank[part]]] is Protected. 
Tag Times in Times[Blank[two], Pattern[two, Blank[part]]] is Protected.
Recursion depth of 1024 exceeded.
Uncaught Throw[TerminatedEvaluation[RecursionLimit], TerminatedEvaluation[RecursionLimit, 1024]] returned to top level.
Tag Times in Times[Blank[Times[-500.0, Plus[Times[0.2, Cos[theta[t]], Power[Derivative[1][phi][t], 
2]], Times[-2, Sin[theta[t]], Derivative[1][theta][t]], Times[0.2, Sin[phi[t]], Derivative[2][phi][t]], Times[2, Cos[theta[t]], Derivative[2][theta][t]]]]], Pattern[two, Blank[part]]] is Protected. 
Tag Times in Times[Blank[two], Pattern[two, Blank[part]]] is Protected.
Recursion depth of 1024 exceeded.
Uncaught Throw[TerminatedEvaluation[RecursionLimit], TerminatedEvaluation[RecursionLimit, 1024]] returned to top level.
3

Never use underscore in the name of variables. Instead of two_part_one write for example twoPartOne. The underscore is used only in patterns.

It is obvious that, when all three functions are identically zero, the system is verified. However, the system in normal form is singular at the origin, that is, the right-hand side contains fractions that become indeterminate for your initial data:

Solve[{(x^\[Prime]\[Prime])[t] == 
    1/2 (-(1/5) Cos[theta[t]] Derivative[1][phi][t]^2 + 
       2 Sin[theta[t]] Derivative[1][theta][t] - 
       1/5 Sin[phi[t]] (phi^\[Prime]\[Prime])[t] - 
       2 Cos[theta[t]] (theta^\[Prime]\[Prime])[t]), (
     theta^\[Prime]\[Prime])[t] == 
    1/3000 (9810 Sin[theta[t]] + 10000 theta[t] + 
       500 (1/5 Cos[2 theta[t]] Derivative[1][phi][t] Derivative[1][
            theta][t] - 
          2 Sin[theta[t]] Derivative[1][theta][t] Derivative[1][x][
            t]) - 500 (1/
           5 Cos[2 theta[t]] Derivative[1][phi][t] Derivative[1][
            theta][t] - 
          2 Sin[theta[t]] Derivative[1][theta][t] Derivative[1][x][
            t] + 1/5 Cos[theta[t]] Sin[theta[t]] (
            phi^\[Prime]\[Prime])[t] + 
          2 Cos[theta[t]] (x^\[Prime]\[Prime])[t])), (
     phi^\[Prime]\[Prime])[t] == 
    1/2100 (-1000 phi[t] - 
       100 Cos[phi[t]] Derivative[1][phi][t] Derivative[1][x][t] - 
       500 (1/5 Cos[theta[t]] Derivative[1][theta][t]^2 + 
          1/5 Cos[theta[t]] Derivative[1][theta][t] Derivative[1][x][
            t] + 1/5 Cos[theta[t]] Sin[theta[t]] (
            theta^\[Prime]\[Prime])[t] + 
          100 Sin[theta[t]] (x^\[Prime]\[Prime])[t]))}, {x''[t], 
   theta''[t], phi''[t]}] // Simplify
% /. t -> 0 /. {x[0] -> 0, x'[0] -> 0, phi[0] -> 0, phi'[0] -> 0, 
  theta[0] -> 0, theta'[0] -> 0 -> Function[0]}
4
Open in Wolfram Cloud ↗