## ode45 unexpected behaviour for initial conditions = 0

Kostas

### Kostas (view profile)

さんによって質問されました 2019 年 10 月 4 日

### Kostas (view profile)

さんによって 編集されました 2019 年 10 月 4 日
Jon

### Jon (view profile)

さんの 回答が採用されました
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)
%% 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.
KMT.

#### 0 件のコメント

サインイン to comment.

## 1 件の回答

### Jon (view profile)

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

Jon

### Jon (view profile)

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
% make overall vector of derivatives
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

### Kostas (view profile)

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

### Kostas (view profile)

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.