Message Boards

WOLFRAM COMMUNITY

228 Views

3 Replies

4 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Circular reference in differential equations

Mishal Mohanlal

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.5m(2ltheta''[t]Cos[theta[t]]-2ltheta'[t]Sin[theta[t]]+2rphi''[t]Sin[phi[t]]+2rphi'[t]^2Cos[theta[t]])-C1x'[t]; two_part_one = -0.5m(2lx''[t]Cos[theta[t]]-2lx'[t]theta'[t]Sin[theta[t]]+2rphi''[t]Cos[theta[t]]Sin[theta[t]]+2rltheta'[t]phi'[t]Cos[2theta[t]]); two_part_two = 0.5m(-2ltheta'[t]x'[t]Sin[theta[t]]+2rltheta'[t]phi'[t]Cos[2theta[t]])+Kbtheta[t]+mglSin[theta[t]]-C2theta'[t]; two = two_part_one+two_part_two; three_part_one = -0.5m(mrx''[t]Sin[theta[t]]+2rx'[t]theta'[t]Cos[theta[t]]+2rltheta''[t]Cos[theta[t]]Sin[theta[t]]+2rltheta'[t]^2Cos[theta[t]]); three_part_two = 0.5m(-2rphi'[t]x'[t]Cos[phi[t]])-Ktphi[t]-C3phi'[t]; three = three_part_one+three_part_two; S = NDSolve[{x''[t] == one/m, theta''[t] == two/(ml^2+2Ib), phi''[t] == three/(mr+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()

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

Sort By:

Mishal Mohanlal

Posted 1 month ago

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.

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

Gianluca Gorni

Gianluca Gorni, University of Udine

Posted 1 month ago

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

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