MATLAB Answers

Error While Solving couple-Differential Equations

62 views (last 30 days)
@ Star Strider
I've tried to solve 5 differential equations. While running the script I'm getting some errors. What shall I modify in the script?
Parameter_TFP.m
function f_Pendulum_structure = Parameter_TFP(t,y, time,U_g, K_s,C_s,m_s,m_1,m_2,m_3,m_4,W,R_eff_1,R_eff_2,R_eff_3,R_eff_4,k_h,mu_1,mu_2,mu_3,mu_4,K_r,d_3,d_2)
% t1 y1 are the variables on which differential equation is formulated
interpolated_Ug = interp1(time,U_g,t);
% friction force calculation
% if condition == presliding
% f_h_1 = K_h *y(1); % pre-sliding condition for concave mass 1
% f_h_2= k_h*(y(3)-y(1)); % pre-sliding condition for concave mass 2
% f_h_3= k_h*(y(5)-y(3)); % pre-sliding condition for concave mass 3
% f_h_4= k_h*(y(7)-y(5)); % pre-sliding condition for concave mass 4
% end
%if condition == sliding
f_h_1= mu_1*W*sign(y(2)); % Sliding condition for concave mass 1
f_h_2= mu_2*W*sign(y(4)-y(2)); % Sliding condition for concave mass 2
f_h_3= mu_3*W*sign(y(6)-y(4)); % Sliding condition for concave mass 3
f_h_4= mu_4*W*sign(y(8)-y(6)); % Sliding condition for concave mass 4
%end
if abs(y(3)-y(1))>d_2
f_r_2 = K_r*(abs(y(3)-y(1))-d_2)*sign(y(3)-y(1));
else
f_r_2 = 0;
end
if abs(y(5)-y(3))>d_3
f_r_3 = K_r*(abs(y(5)-y(3))-d_3)*sign(y(5)-y(3));
else
f_r_3 = 0;
end
f_Pendulum_structure = [y(2);(-interpolated_Ug) + (f_h_2/m_1)-(f_h_1/m_1)-(W/(m_1*R_eff_1))*y(1)+((W/(R_eff_2*m_1))*(y(3)-y(1)))+(f_r_2/m_1);...
y(4);(-interpolated_Ug) + (f_h_3/m_2)-(f_h_2/m_2)-((W/(m_2*R_eff_2))*(y(3)-y(1)))+((W/(R_eff_3*m_2))*(y(5)-y(3)))+(f_r_3/m_2)-(f_r_2/m_2);...
y(6);(-interpolated_Ug) + (f_h_4/m_3)-(f_h_3/m_3)-((W/(m_3*R_eff_3))*(y(5)-y(3)))+((W/(R_eff_4*m_3))*(y(7)-y(5)))-(f_r_3/m_3);...
y(8);(-interpolated_Ug) - (f_h_4/m_4)-((W/(m_4*R_eff_4))*(y(7)-y(5)))+(C_s*(y(10)-y(8))/m_4)+(K_s*(y(9)-y(7))/m_4);...
y(10);(-interpolated_Ug) -(C_s*(y(10)-y(8))/m_s)-(K_s*(y(9)-y(7))/m_s)];
end
main file
clear all;
close all;
clc
%--------------------------------------------------------------------------
tic
%--------------------------------------------------------------------------
%% Ground Acceleration Data------------------------------------------------
load EQ_data_Peak.txt;
skip = 1;
time= EQ_data_Peak(1:skip:end,1);
U_g= 9.81*EQ_data_Peak(1:skip:end,2);
TIME_STEP = time(2) - time(1);
%%-------------------------------------------------------------------------
%%%%%%%%%%% Friction Pendulum with Structure as SDOF system
% model parameters
freq_natural = 1.2; % in Hz
T_natural = 1/freq_natural; % in sec
m_s = 1000*1223.2; % modal mass of the structure in kg
K_s = ((2*pi*freq_natural).^2)*m_s; %%%% stiffness of structure in N/m
eta = 0.01; % damping ratio
C_s = 2*eta* sqrt(K_s*m_s); % damping co-efficient of structure
omega = sqrt(K_s/m_s); % eigen-frequency (rad/s)
omega_d = omega.*sqrt(1-eta.^2); % damped eigen frequency (rad/s)
m_1 = 465; % in kg
m_2 = 465; % in kg
m_3 = 258; % in kg
m_4 = 853; % in kg
W = 12000*1000 ; % in N
R_eff_1 =1.522 ;% in m
R_eff_2 =0.190 ;% in m
R_eff_3 =0.190 ;% in m
R_eff_4 =1.522 ;% in m
k_h = (W/R_eff_4)+100; % Pre-sliding stiffness co-efficient
mu_1= 0.02;
mu_2 =0.0175;
mu_3 =0.0175;
mu_4 =0.02;
K_r = (W/R_eff_4)+100;
d_3 =0.07; % in m
d_2 =0.07; % in m
%%%%%%% initial condition provided
initial_condition = [0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0]; %values for start of Triple Pendulum Bearing along with SDOF Structure
%% Solving for 4 Concave Bearings Plate Along with SDOF Structure
[t,y] = ode15s(@(t,y) Parameter_TFP(t,y,time,U_g, K_s,C_s,m_s,m_1,m_2,m_3,m_4,W,R_eff_1,R_eff_2,R_eff_3,R_eff_4,k_h,mu_1,mu_2,mu_3,mu_4,K_r,d_3,d_2), time,initial_condition);
toc
error:
Error using vertcat
Dimensions of matrices being concatenated are not consistent.
Error in Parameter_TFP (line 29)
f_Pendulum_structure = [y(2);(-interpolated_Ug) + (f_h_2/m_1)-(f_h_1/m_1)-(W/(m_1*R_eff_1))*y(1)+((W/(R_eff_2*m_1))*(y(3)-y(1)))+(f_r_2/m_1);...
Error in Structure_TFP_Bearing>@(t,y)Parameter_TFP(t,y,time,U_g,K_s,C_s,m_s,m_1,m_2,m_3,m_4,W,R_eff_1,R_eff_2,R_eff_3,R_eff_4,k_h,mu_1,mu_2,mu_3,mu_4,K_r,d_3,d_2)
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 150)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in Structure_TFP_Bearing (line 46)
[t,y] = ode15s(@(t,y) Parameter_TFP(t,y,time,U_g, K_s,C_s,m_s,m_1,m_2,m_3,m_4,W,R_eff_1,R_eff_2,R_eff_3,R_eff_4,k_h,mu_1,mu_2,mu_3,mu_4,K_r,d_3,d_2),
time,initial_condition);

Accepted Answer

Walter Roberson
Walter Roberson on 3 Mar 2021
y(10);(-interpolated_Ug) -(C_s*(y(10)-y(8))/m_s)-(K_s*(y(9)-y(7))/m_s)];
You have the interpolation variable, space minus-sign expression
The previous line you have the interpolation variable space minus-sign space expression. The second space makes a difference.
Consider the difference between
[1 - 2]
[1 -2]
The first one is subtraction, one minus two. The second one is a list with two elements, one and negative two.
  1 Comment
Sumit Saha
Sumit Saha on 3 Mar 2021
@Walter RobersonThe code is running for 2 hr. Why it's taking so much time ? Is there any option or way to rewrite the code in a different manner?

Sign in to comment.

More Answers (2)

William Rose
William Rose on 3 Mar 2021
In function Parameter_TFP(), 10 rows are combined to make vector f_Pendulum_structure. Matlab thinks the rows being combined do not all have the same length, which causes a vertcat() error. I don;t know why Matlab thinks this, since they all appear to be scalars, but I fixed this error by putting a squeeze() around each row in the lines that assign a value to f_Pendulum_structure. I have found that helpful before.
Also, the call to ode15s() passes the parameter time which is (6000x1). This should be tspan=[tstart, tend], which is (1x2). I fixed that by adding tspan and passing it, instead of time, to ode15s(). (time is still passed to Parameter_TFP.)
With these fixes, I no longer get an error in Parameter_TFP, but I get a new error:
>> main_wr
Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the
step size below the smallest value allowed (7.905050e-323) at time t.
> In ode15s (line 668)
In main_wr (line 48)
Elapsed time is 0.442605 seconds.
I changed ode15s() to ode45(). That eliminated the warning error above, but it seems to be very slow, because it did not finish after several minutes. I reduced the end time from 29.995 sec to 0.1 sec, and still, main_wr did not complete in the minute that I waited. But it didn't give any errors either. This suggests to me that you need to examine your differential equations. I notice that every time I control-C out of main_wr (because it is taking forever), Matlab says it was interrrupted at line 3 of the modified Parameter_TFP_wr() function. This is the line that calls interp1(). Does that mean line 3 is slow? Could line 3 be eliminated or changed? It is not obvious that it could be, since you want U_g at the current time, what ever that is. The modified main() and the modified Parameter_TFP() are attached. main_wr() does not generate errors, but I have not seen it finish.
  13 Comments
William Rose
William Rose on 6 Mar 2021
@Walter Roberson, that is inersting about what's in the finance toolboxes. Your knowledge of Matlab and its many toolboxes is amazing.
One of my children took, and passed, an actuary exam that included stochastic differential equations: differential equations that describe how probability distributions evolve over time. I'm glad I wasn't the one taking that exam!

Sign in to comment.


William Rose
William Rose on 4 Mar 2021
I ran your code, which includes non-smooth functions that are bad for ODE solvers, as @Walter Roberson explained: four statemnts call sign(), for f_h_1 to f_h_4, and there are two if statements, which also include calls to abs() and sign(), which are used to compute f_r_2 and f_r_3.
I added two smooth functions to replace the non-smooth ones: smoothSign(x,xw) and smoothHeavi(x,xw). smoothSign is a substitute for sign(x). smoothHeavi(x,xw) is a smooth version of the Heaviside function, which exists in the Symbolic Math toolbox but not in regular Matlab. By the way, if you specify a very narrow transition width xw, then x/xw or -x/xw can exceed 1000, and then exp(x/xw) returns NaN. My smoothSign(x,xw) and smoothHeavi(x,xw) protect against that.
When I integrated with ode45() from t=0 to 0.05, without the smooth functions, the execution time was 268 sec and the integrator took 24.5 million steps.
When I replaced the non-smooth statements with calls to smoothSign() and smoothHeavi(), the integraiton did not get faster and the number of steps did not shrink, as I had hoped, at least not whn xw=1e-8 or 1e-6. But when I increase the transition width for smoothSign to 1e-4 (i.e. smoother function), the executon time dropped by a factor of 200 and the number of steps decreased by a factor of 40. That is progress - if the results are acceptable. You will have to experiment with the values of yw and xw and examine results to see if they seem plausible. Maybe ode45() is not the best solver for this problem. I am attaching main_wr2.m and Parameter_TFP_wr2.m. The smooth functions are included in Parameter_TFP_wr2.m.
  1 Comment
Sumit Saha
Sumit Saha on 4 Mar 2021
Thanks a lot...May be ode 23s can work...I'll check that and let you know regarding this issue.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!