I have some coupled nonlinear ordinary differential equations. Three equations have three 2nd order derivative coupled together. Can someone tell me how I can get the response

I have some coupled nonlinear ordinary differential equations. Three equations have three 2nd order derivative coupled together. Can someone tell me how I can get the response, numerically? i.e t vs x , t vs theta dot plots.

Thanks in advance

[James Tursa]

Do you have initial conditions for all the variables and first derivatives?

[prabin kumar meher]

all initial condition are zero

Answer 1

[Torsten]

Solve the third equation for theta_dotdot and insert this expression in equations (1) and (2).

Then use the usual substitutions

z1 = x
z2 = xdot
z3 = y
z4 = ydot
z5 = theta
z6 = thetadot

[prabin kumar meher]

Torsten thanks for your response , if there is slight change in the third equation, then what should i do.

[James Tursa]

This is not a slight change, it is a huge change. By including x_dotdot and y_dotdot in that third equation you can no longer use Torsten's advice directly because these terms would show up inside the f( ) function argument, preventing you from solving all the equations for the highest order derivatives. What does f( ) look like?

[prabin kumar meher]

so @James Tursa please tell me how shold i proceed to solve via ODE45.

[Torsten]

@James Tursa

In my opinion, it's a constant named "f". But the OP will clarify, I guess.

@prabin kumar meher

In this case, solve the first equation for x_dotdot and the second for y_dotdot.

Insert x_dotdot and y_dotdot in the last equation and solve for theta_dotdot. Then insert the expression you obtained for theta_dotdot in the first and second equation and proceed as before.

[prabin kumar meher]

yes f is a constant term

[James Tursa]

OK, then you can use Torsten's advice.

Answer 2

[James Tursa]

Another related approach is to isolate all the 2nd derivatives on the LHS and put everything else on the RHS. Then pick off the coefficients of the 2nd order derivatives to form a matrix problem that you solve for the 2nd derivatives. E.g., the matrix problem would look something like this based on your comment edits:

A*z = b

where

A = [     1                0         -f*sin(theta);
          0                1          f*sin(theta);
    -m*e*sin(theta)  m*e*cos(theta)  (m*e^2+Ip)/B^2]

and

z = [x_dotdot; y_dotdot; theta_dotdot]

and

b = 3x1 vector with all the other stuff from RHS

Then solve for z to get the 2nd derivatives.

[prabin kumar meher]

Warning: Failure at t=1.815139e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed

(5.684342e-14) at time t.

above warning is arise. what should i do now.

[Torsten]

For that we are able to add something useful, you must provide your code.

[prabin kumar meher]

clc;
close all;
clear all;

% T=0;0.01;100;
% IC=[0 0.1 0 0.1 0 0.1];
t=0:0.1:90;
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,x]=ode45(@system_prabin,t,[0 0 0 0 0 0 ],options);
plot(t,x(:,5))

function dx = system_prabin(t,x)
e=10*5*10^-7;
m=0.54;
C=0.0458;
C0=0.5*10^-3;
K=4.865*10^-7;
mum=0.1;
Rd=0.0002;
Rm=40;
Ip=0.1;
Kp=2440;
Kd=0.1178;
I0=0.8;
Vs=1;
a=0.3927;

f=e/C0;
p=Kp*C0/I0;
B=sqrt((m*C0^3)/(K*I0^2));
d=Kd*C0/(B*I0);
c1=B*C/m;
% P=(mum*Vs)/(Rm*((m*e^2)+Ip));
% S=((mum^2)/(Rm*((m*e^2)+Ip)))+((Rd)/((m*e^2)+Ip));
mu= c1+(8*d*(cos(a)+cos(3*a)));
w= sqrt(8*p*(cos(a)+cos(3*a))-16);
a1=-24-4*cos(4*a)-4*cos(12*a)+6*p*(3*cos(a)+4*cos(3*a)+cos(9*a))-4*p^2*(cos(2*a)+cos(6*a)+2);
a2=-24+12*cos(4*a)+12*cos(12*a)+18*p*(cos(a)-cos(9*a))+4*p^2*(cos(2*a)+cos(6*a)-2);
a3=-8*p*d*(cos(2*a)+cos(6*a)+2)+6*d*(3*cos(a)+4*cos(3*a)+cos(9*a));
a4=6*d*(cos(a)-cos(9*a));
a5=4*d^2*(cos(2*a)+cos(6*a)-2);
a6=-4*d^2*(cos(2*a)+cos(6*a)+2);
a7=12*d*(cos(a)-cos(9*a))+8*p*d*(cos(2*a)+cos(6*a)-2);
dx= zeros(6,1);
dx(1)=x(2);
dx(2)=(f*((x(6)^2*cos(x(5)))+(((1/(((m*e^2)+Ip)-(f*m*e*sin(x(5))*sin(x(5)))-(f*m*e*cos(x(5))*sin(x(5)))))*(((mum*Vs*B^2)/Rm)-(((mum^2/Rm)+Rd)*B*x(6))+(f*m*e*cos(x(5))*(sin(x(5))-cos(x(5)))*x(6)^2)-(m*e*sin(x(5))*((mu*x(2))+(w^2*x(1))+(a1*x(1)^3)+(a2*x(1)*x(3)^2)+(a3*x(1)^2*x(2))+(a4*x(2)*x(3)^2)+(a5*x(1)*x(4)^2)+(a6*x(1)*x(2)^2)+(a7*x(1)*x(3)*x(4))))+(m*e*cos(x(5))*((mu*x(4)-(w^2*x(3))+(a1*x(3)^3)+(a2*x(3)*x(1)^2)+(a3*x(3)^2*x(4))+(a4*x(4)*x(1)^2)+(a5*x(3)*x(2)^2)+(a6*x(3)*x(4)^2)+(a7*x(1)*x(3)*x(2)))))))*sin(x(5)))))-(mu*x(2))-(w^2*x(1))-(a1*x(1)^3)-(a2*x(1)*x(3)^2)-(a3*x(1)^2*x(2))-(a4*x(2)*x(3)^2)-(a5*x(1)*x(4)^2)-(a6*x(1)*x(2)^2)-(a7*x(1)*x(3)*x(4));
dx(3)=x(4);
dx(4)=(f*((x(6)^2*cos(x(5)))-(((1/(((m*e^2)+Ip)-(f*m*e*sin(x(5))*sin(x(5)))-(f*m*e*cos(x(5))*sin(x(5)))))*(((mum*Vs*B^2)/Rm)-(((mum^2/Rm)+Rd)*B*x(6))+(f*m*e*cos(x(5))*(sin(x(5))-cos(x(5)))*x(6)^2)-(m*e*sin(x(5))*((mu*x(2))+(w^2*x(1))+(a1*x(1)^3)+(a2*x(1)*x(3)^2)+(a3*x(1)^2*x(2))+(a4*x(2)*x(3)^2)+(a5*x(1)*x(4)^2)+(a6*x(1)*x(2)^2)+(a7*x(1)*x(3)*x(4))))+(m*e*cos(x(5))*((mu*x(4)-(w^2*x(3))+(a1*x(3)^3)+(a2*x(3)*x(1)^2)+(a3*x(3)^2*x(4))+(a4*x(4)*x(1)^2)+(a5*x(3)*x(2)^2)+(a6*x(3)*x(4)^2)+(a7*x(1)*x(3)*x(2)))))))*sin(x(5)))))-(mu*x(4))-(w^2*x(3))-(a1*x(3)^3)-(a2*x(3)*x(1)^2)-(a3*x(3)^2*x(4))-(a4*x(4)*x(1)^2)-(a5*x(3)*x(2)^2)-(a6*x(3)*x(4)^2)-(a7*x(1)*x(3)*x(2));
dx(5)=x(6);
dx(6)=(1/(((m*e^2)+Ip)-(f*m*e*sin(x(5))*sin(x(5)))-(f*m*e*cos(x(5))*sin(x(5)))))*(((mum*Vs*B^2)/Rm)-(((mum^2/Rm)+Rd)*B*x(6))+(f*m*e*cos(x(5))*(sin(x(5))-cos(x(5)))*x(6)^2)-(m*e*sin(x(5))*((mu*x(2))+(w^2*x(1))+(a1*x(1)^3)+(a2*x(1)*x(3)^2)+(a3*x(1)^2*x(2))+(a4*x(2)*x(3)^2)+(a5*x(1)*x(4)^2)+(a6*x(1)*x(2)^2)+(a7*x(1)*x(3)*x(4))))+(m*e*cos(x(5))*((mu*x(4)-(w^2*x(3))+(a1*x(3)^3)+(a2*x(3)*x(1)^2)+(a3*x(3)^2*x(4))+(a4*x(4)*x(1)^2)+(a5*x(3)*x(2)^2)+(a6*x(3)*x(4)^2)+(a7*x(1)*x(3)*x(2))))));

end

[Torsten]

Works under R2022 b (see above).

[prabin kumar meher]

yes@Torsten it is working up to 20 sec ,but as i increase the time range to 1000, it is not working.

[Torsten]

It works up to t = 90 approximately (see above).

Then all your variables seem to blow up.

Either your equations are wrong, there is a singularity around t = 90 or the integrator is not able to solve your problem.

For all three problems, I cannot be of help.

[prabin kumar meher]

Thank you @Torsten for yor help