MATLAB Answers

Kostas
0

ode45 unexpected behaviour for initial conditions = 0

Kostas
さんによって質問されました 2019 年 10 月 4 日
最新アクティビティ Kostas
さんによって 編集されました 2019 年 10 月 4 日
Hi all,
I have a code where I am using it to solve the following differential equation:
In specific, the above is a 4x4 system of equations where J is the polar moment of inertia (a diagonal matrix), K is a stiffness matrix, and T_g(th) is the excitation torque on the system.
What I want to do is simply solve the system for some initial conditions that I set myself. However, just as a sanity check, I solve the equation with all initial conditions set to zero. Naturally, this should give a zero response for the equation's solution, however it doesn't which I find quite odd. Hence with that result, I cannot trust that I will have a correct answer with any initial condition.
So any ideas of what might be going wrong?
Attached you can find the .mat file that I draw all my inputs from and my ode45 code is below:
clear
clc
% Input data
% Input Engine Data
omega = 94 * pi / 30 ; % Rotational Speed as initial condition
% Degrees of Freedom
do = [1 4 7 11 18] ; % DOF Split vector
N = length(do) - 1 ; % Number of DOFs
% ODE Solution
tmax = 50 ; % Solution time (s)
% Load ODE Matrices
load('Struct.mat')
%% ODE Solution
% ODE Initial Conditions
% Time
tspan = [0 tmax];
% Initial conditions
x0 = [zeros(1,N) repmat(omega, 1, N)] ;
% ODE options
options = odeset('Mass', [eye(N) zeros(N) ; zeros(N) J]) ;
% Solution
[t, xSol] = ode45(@(t, th) odefcn(t, th, thx, K, Tg, N) , tspan, x0, options) ;
%%
figure
subplot(2, 1, 1)
plot(t, xSol(:,1:N))
title('Displacement')
subplot(2, 1, 2)
plot(t, xSol(:,N+1:end))
title('Velocity')
%% ODE
function dxdt = odefcn(~, th, thx, K, Tg, N)
% Indices
m = N + 1 ;
n = N * 2 ;
% Preallocate
dxdt = zeros(n, 1) ;
% Torques
T = diag(interp1(thx, Tg', wrapToPi(th(1:N)))) ; % Gas Torque
% ODE Equation
dxdt(1:N) = th(m:n) ;
dxdt(m:n) = - K * th(1:N) + T;
end
P.S.
I have tried different solvers and I am getting the same unexpected result.
Thanks for your help in advance,
KMT.

  0 件のコメント

サインイン to comment.

1 件の回答

Jon
回答者: Jon
2019 年 10 月 4 日
 採用された回答

I think you have a basic conceptual error in your ode definition.
You should have an overall state vector x = [x1;x2] whre x1 = theta, x2 = d(theta)/dt so your ode's look like
dx1/dt = x2
dx2/dt = -K*theta + T(x)
you have instead dx1/dt = theta, which is not correct

  4 件のコメント

Jon
2019 年 10 月 4 日
Yes, I see that it looks like your equations are equivalent to mine. I didn't look carefully enough at your indices. For the future, the code might be a little more obvious if you did something like
function dxdt = odefcn(~, th, thx, K, Tg, N)
% Torques
T = diag(interp1(thx, Tg', wrapToPi(th(1:N)))) ; % Gas Torque
% ODE Equation
% assign local variables for readability
theta = th(1:N);
omega = th(N+1:end);
% compute derivatives
dthetadt = omega;
domegadt = -K*theta + T;
% make overall vector of derivatives
dxdt = [dthetadt;domegadt];
When you do your test setting all of the initial conditions to zero, do you somehow also force the torque to zero. Otherwise, if I understand your system correctly, wouldn't a non-zero torque cause the system to accelerate even if the initial conditions were zero?
Kostas
2019 年 10 月 4 日
''...I understand your system correctly, wouldn't a non-zero torque cause the system to accelerate even if the initial conditions were zero?''
This is a valid point, and I have been thinking of that as well, so here is my logic: The torque is strictly a function of theta and not a function of time (which is the independent variable). In addition, the torque function is such that T(0)=0. In which case I figure that if theta begins as zero, there will be no reason for the excitation torque to become anything different than zero as well. Hence the response should be zero as well. By the same token, if the excitation torque was time dependent, then yes as you said the response would be non-zero as well.
Kostas
2019 年 10 月 4 日
Indeed you are correct however! I just had a look and as the function is approximated near zero, it is not exactly equal to zero when it should. Given how small the aplitudes of the oscillation are I would say that this is what could be causing this

サインイン to comment.



Translated by