MATLAB Answers

Non-linear coupled ODE system of equations

조회 수: 34(최근 30일)
Panagiota Atzampou
Panagiota Atzampou 1 Feb 2021
답변: Jon 1 Feb 2021
Greetings,
I want to solve the following set of coupled ODEs:
I went ahead and used the following assumptions: [ y(1)=θ(t), y(2)=d(θ(t))/dt ,y(3)=φ(t) and y(4)=d(φ(t))/dt ] in order to later on employ the ode45 command to find the solution. Eventhough the script I made runs with no errors or issues, the output matrix y is filled with NaN (except the first row of course, which corresponds to the initial conditions given).
Could you please help me pin point the error in my code, or alternatively propose another method perhaps more suitable to solve these equations?
Thank you in advance!
Here's the script:
% Definition of Parameters
g=9.81; % m/s2
l= input('Strings length: ');
omega=(g/l);
Ttot=input('Total Simulation Time (sec): ');
% Initial Conditions
theta_0 = input('Theta_0 (in form of n*pi): ');
theta_t0 = input('Theta_t0: ');
phi_0 = input('Phi_0 (in form of n*pi): ');
phi_t0 = input('Phi_t0: ');
y0 = [theta_0, theta_t0, phi_0, phi_t0];
% Solution
f=@(t,y) Pendulum3D(t,y,omega);
[t,y]=ode45(f,[0 Ttot],y0);
% Function
function dy=Pendulum3D(t,y,omega)
dy=[y(2);(sin(y(1))*cos(y(1))*(y(4))^2)-omega*sin(y(1)); y(4); (-2*y(4)*y(2)*cot(y(1)))];
end
  댓글 수: 1
Bjorn Gustavsson
Bjorn Gustavsson 1 Feb 2021
Check what happens when cot(theta) goes to infinity, that is the first worrisome point in your code.

댓글을 달려면 로그인하십시오.

채택된 답변

Jon
Jon 1 Feb 2021
As Bjorn has suggested, you must guard against values of theta where cot(theta) goes to +/- infinity, e.g. 0,pi
Your code seems to run ok for the initial conditions I've edited in below (I commented out the prompts for user input and just supplied values so the results are reproducible)
% Definition of Parameters
g=9.81; % m/s2
l= 1; %input('Strings length: ');
omega=(g/l);
Ttot=10;% input('Total Simulation Time (sec): ');
% Initial Conditions
theta_0 = pi/2; %input('Theta_0 (in form of n*pi): ');
theta_t0 = 0; % input('Theta_t0: ');
phi_0 = pi/4; %input('Phi_0 (in form of n*pi): ');
phi_t0 = 1; %input('Phi_t0: ');
y0 = [theta_0, theta_t0, phi_0, phi_t0];
% Solution
f=@(t,y) Pendulum3D(t,y,omega);
[t,y]=ode45(f,[0 Ttot],y0);
% plot results
plot(t,y)
I would suggest making your equations more obvious as follows, but it seems that yours are equivalent, just hard to read
function ydot = Pendulum3D(~,y,omega)
% use actual variable names to make code more understandable, help avoid
% errors in formulas
theta = y(1);
thetaDot = y(2);
phi = y(3);
phiDot = y(4);
ydot(1) = thetaDot;
ydot(2) = sin(theta)*cos(theta)*phiDot^2 - omega*sin(theta);
ydot(3) = y(4);
ydot(4) = -2*phiDot*thetaDot*cot(theta);
% make sure it returns as a column
ydot = ydot(:);
end

추가 답변(0개)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by