how to solve the problem of ode45 ?

Dear,

I have been working on a simple model. I used the ode45 solver, and one of my outputs doesn't make sense, I was wondering if is that due to the ode45 or a convergence-related issue. Can you please have a look at the very simple code attached here?

the question is that in the code, once I use C_n2=15*10^7, the result ("ALLACC") was supposed to be smaller than those of the C_n2=10^7 result ("ALLACC").

Is that due to ode45? or is it due to something else?

Can you please help me with that?

That would be much appreciated if you could help me.

for running the code, I attached the load file (INPUT.txt) too.

thanks again for your tuime.

Regards,

funcfile:

function [xdot]=funcfile(t,x, C_n2)    

xdot=zeros(4,1);   

global acc tspan m1 m2  k1 k2 c1 c2 

z_t=interp1(tspan,acc,t,'spline');

%% [M] 
m1=4.9814*10^6;
m2=3.4382*10^6; 
%% [K]
k1=1.042585220149801*10^8; 
k2=2.232017093496502*10^9; 

%% [C]
c1=(1.065129353782337*10^7)+5.794706239147969*10^-4; 
c2=1.443681417713320*10^7; 

NOCo2=x(2)*x(4);

if NOCo2<0

    NOD2=C_n2*x(4);

else

   NOD2=0*x(4);  

end


M=[m1 0;
    0 m2];

C=[c1+c2 -c2;
   -c2 c2];

K=[k1+k2 -k2;
   -k2 k2];

Acc_m=[-m1 -m2]';

F=Acc_m.*z_t;

NOF=[0 NOD2]';


xdot(1:2)=x(3:end);
xdot(3:end)=inv(M)*(-C*x(3:end)-K*x(1:2)-NOF+F);

end 

Mainfile:

clear all
clc

global acc tspan  m2  k2  c2  

dt=0.01;

load('INPUT.txt');

acc= (INPUT(:,2))';

t = (0:(length(acc)-1))*(dt);
t=(t');
u=[acc]; 
tspan=t';

z_t=interp1(tspan,acc,t,'spline');

x0=zeros(1,4);

C_n2=10^7;  % first 

% C_n2=15*10^7;  % second
 
[t,x]=ode45(@funcfile,tspan,x0,[],C_n2 );
 
NOCo2 = x(:,2).*x(:,4);
NOD2 = C_n2*x(:,4).*(NOCo2 < 0) + 0*x(:,4).*(NOCo2 >= 0);

ACC2=(-m2.*z_t- NOD2-((-c2).*x(:,3)+c2.*x(:,4))- ...
     ((-k2).*x(:,1)+k2.*x(:,2)))./m2;

ALLACC=[ max(abs(ACC2)) ];

[Star Strider]

I cannot understand what you want to do from your code. If you intend to estimate the parameters of whatever ‘funcfile’ describes to your data, there are ways to do that (however that may not be straightforward with this file, even assuming it is possible). Also, since the parameters vary over about 12 orders of magnitude, the system is considered to be ‘stiff’ and ode15s may be a better choice to integrate the differential equatiions than ode45.

It would also be helpful to share the symbolic form of the equations that are coded in ‘funcfile’ as well as a description of what you want to do and the parameters you want to estimate (if that is what you want to do).

A1 = readmatrix('INPUT.txt');

figure
plot(A1(:,1), A1(:,2))
grid
xlabel('Column #1')
ylabel('Column #2')
title('INPUT.txt Data')

.

[ADNAN]

many thanks for your reply @Star Strider, I named some parameters in the file and updated that

in the equation of motion of the system, there is a force named "NOD2"(it is a viscous damper force), and "C_n2" is the damping coefficient.

As you know ode45 can give us the displacement (x) and velocity (xdot) of a 2-degree-of-freedom system. After the simulation is completed, I want to calculate the acceleration of the system named as "ACC2" in Mainfile . Finally, I want to check maximum of that (ALLACC=[ max(abs(ACC2)) ]).

I did the simulation under two damping coefficients ("C_n2=10^7" and "C_n2=15*10^7"), for "C_n2=15*10^7", I was supposed to get smaller acceleration (ALLACC=[ max(abs(ACC2)) ]) than once I use "C_n2=10^7", yet I did not get. I got bigger "ALLACC" for "C_n2=15*10^7".

I have checked every detail but there was not any problem, I just thought that it could be due to the solver.

Can you please have a look is that due to matlab solution ? or ?

many thanks for your time again.

Kind regards,

Mainfile:

clear all
clc
global acc tspan  m2  k2  c2  

dt=0.01; %  time step 
load('INPUT.txt'); % accelaration input 

acc= (INPUT(:,2))';    % accelaration input 
t = (0:(length(acc)-1))*(dt); 
t=(t');
u=[acc];  % accelaration input 
tspan=t'; 
z_t=interp1(tspan,acc,t,'spline'); % new acceleration after interpolations

x0=zeros(1,4);

C_n2=10^7;  % first damping coefficient 
% C_n2=15*10^7;  % second damping coefficient 
 
[t,x]=ode45(@funcfile,tspan,x0,[],C_n2 );
 
% calculation of damping force after the simulation

NOCo2 = x(:,2).*x(:,4); % relation of force
NOD2 = C_n2*x(:,4).*(NOCo2 < 0) + 0*x(:,4).*(NOCo2 >= 0); % calculation of damper force

% calculation of accelaration of the system after simulation completed

ACC2=(-m2.*z_t- NOD2-((-c2).*x(:,3)+c2.*x(:,4))- ...
     ((-k2).*x(:,1)+k2.*x(:,2)))./m2;

% calcualtion of maximum absoulte accelation of the 2-degree-of-system

ALLACC=[ max(abs(ACC2)) ];

funcfile:

function [xdot]=funcfile(t,x, C_n2)    
xdot=zeros(4,1);   
global acc tspan m1 m2  k1 k2 c1 c2 
z_t=interp1(tspan,acc,t,'spline');
%% [M] mass 
m1=4.9814*10^6;
m2=3.4382*10^6; 
%% [K] stiffness
k1=1.042585220149801*10^8; 
k2=2.232017093496502*10^9; 
%% [C]  damping 
c1=(1.065129353782337*10^7)+5.794706239147969*10^-4; 
c2=1.443681417713320*10^7; 

% damping force calculation 

NOCo2=x(2)*x(4);
if NOCo2<0
    NOD2=C_n2*x(4);
else
   NOD2=0*x(4);  
end


% mass matrix 
M=[m1 0;
    0 m2];

%damping matrix
C=[c1+c2 -c2;
   -c2 c2];

%stiffness matrix 

K=[k1+k2 -k2;
   -k2 k2];
% mass
Acc_m=[-m1 -m2]';

%force
F=Acc_m.*z_t;

NOF=[0 NOD2]';

xdot(1:2)=x(3:end);
xdot(3:end)=inv(M)*(-C*x(3:end)-K*x(1:2)-NOF+F);
end 

[Star Strider]

I am still lost.

The solver is not likely to be the problem, although for differential equations with several orders-of-magnitude in their coefficients, a ‘stiff’ solver is the best choice.

If you want to estimate the parameters of the differential equations for your data, see Monod kinetics and curve fitting - MATLAB Answers - MATLAB Central or Parameter estimation for a system of differential equations with multiple time spans - MATLAB Answers - MATLAB Central or Reaction kinetics parameter determination - MATLAB Answers - MATLAB Central or any number of others posted in the last decade.

Also, please do not use global variables, since that is not considered to be good programming practice for a number of reasons, not the least of which is that functions using them can change them, causing problems that can be next to impossible to detect and solve. Pass them as extra parameters instead, as described in Passing Extra Parameters.

Consider using the mldivide, \ function here, instead of inv, changing:

xdot(3:end) = inv(M)*(-C*x(3:end)-K*x(1:2)-NOF+F);

to:

xdot(3:end) = M\(-C*x(3:end)-K*x(1:2)-NOF+F);

.

[Walter Roberson]

if NOCo2<0

You need to use an event function and terminate integration. ode*() functions will produce wrong results for discontinuities.

[ADNAN]

thanks for your reply @Walter Roberson, can you please rewrite that part for me ? I could not figure it out. Thanks again.

[Walter Roberson]

See ballode example.