ode45 unexpected behaviour for initial conditions = 0
3 次查看(过去 30 天)
显示 更早的评论
Hi all,
I have a code where I am using it to solve the following differential equation:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/241153/image.png)
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 个评论
采纳的回答
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 个评论
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!