Message Boards Message Boards

Circular reference in differential equations

Posted 1 month ago

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()
POSTED BY: Mishal Mohanlal
3 Replies

POSTED BY: Michael Rogers

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]}
POSTED BY: Gianluca Gorni

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.
POSTED BY: Mishal Mohanlal
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