Error While Solving couple-Differential Equations

@ 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);
 

Answer 1

[Walter Roberson]

           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.

[Sumit Saha]

@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?

Answer 2

[William Rose]

@Sumit Saha

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.

[Sumit Saha]

Thanks @William Rose . Actually through U_g I'm just calling the data at each time step. I've also run the program using ode 15s and got the same error. Now I'm running with the ode 23s. @Walter Roberson what is you suggestion on this issue ?

[Sumit Saha]

@William Rose did you excute the scripts ?

[William Rose]

@Sumit Saha

I executed the scripts you originally posted and got the same errors you reported. I executed the scripts I posted for you, in the posted form, and also some intermediate versions. For example, I tried ode15s() and ode45(), and I posted it with ode45().

[Sumit Saha]

@William Rose do you have any idea about the alternate solving procedure?

[William Rose]

@Sumit Saha

I applied the fix to the spacing issue, which @Walter Roberson astutely identified.

I changed the U_g_interp from interpolation on U_g to a simple lookup from U_g, rounding to the nearest time point in U_g, in an attempt to save time. I plotted the interpolated and looked-up values and the difference is small.

I ran the program successfully to completetion with ode45(), with an end time of 0.1 sec. It took 703 seconds and it produced 59 million time points, i.e. a point every 17 nanoseconds, which is a crazy high rate. No wonder it is taking so long. I did

>> histogram(t) 

to view a plot of the distribution of time points. They were reasonably distributed across the entire 0.1 s interval, somewhat more densely packed in the center region of the integration period. I had wondered if the ingtegrator was taking many fine steps at the start, to get going, and then takign bigger steps later, but the histogram plot showed that was not the case.

I added a progress indicator: seconds completed this far in the integration, and I reran the program, with an end time of 0.03 seconds. It took 115 seconds and produced 11 million time points. See output below:

>> main_wr2
0.005 0.010 0.015 0.020 0.025 
Elapsed time is 135.274376 seconds.
>> disp(size(y))
    11465657          10

File main_wr2.m calls Parameter_TFP_wr2.m. Both files are attached.

I think the integrator is taking many small steps because you have two if statements in Parameter_TFP. , for f_r_2 and f_r_3. I assume thse have to do with static and sligin friction but I am not sure. I have had this problem myself when simulating flow thorugh one-way vaves in the heart. A one-way valve has zero flow when the prssure gradient is negative and a a linear pressure-flow relation when the gradient is positive. I simulated this with an if statement, and the integrator ground to a halt. It basically chatters, taking many tiny time steps when the pressure is near the forward/reverse bias point. I cured it by replacing the if statement with a smooth function that turns on gradually. I think you should do likewise. I also notice that one arm of each if statement includes abs() and sign(), which are als discontinuous functions. We may need to adjust those also. Please explain how fr2 is supposed to depend on y1, y3, and d2, and why. Provide plots if that will help. This explanaton could help design a smooth function to replace the if statements and maybe the abs() and sign() statements.

[Walter Roberson]

The ode*() routines all assume that the functions implemented are C2 (continuous second derivative): this is part of the mathematics of the theory used.

This implies that you should not have any if in the functions unless the if always takes the same branch for any one call to the ode*() routines, or the code has been carefully designed so that the boundary points have continuous second derivatives. In situations where that does not hold, you should be using event functions to terminate the integration and then call the function again to restart it on the other side of the boundary. The ode*() routines will often not notice if the function is only C1, but the results you get may be wrong. If the functions do not have a continuous first derivative, then the ode*() routines are more likely to notice -- but it depends on how large the discontinuity is. If the size of the discontinuity is less than a quadratic projection of the current error margin, then it might not notice... but again the results you get may be wrong.

When I say "might not notice" in this context, I mean that it might not explicitly stop saying that it could not meet integration tolerances. But it might take a long time to do the integration.

[William Rose]

@Sumit Saha,@Walter Roberson

Here is a function that does the if statement thing

if x>0, y1=x; else, y=0; end

in a smooth way:

xw=0.2;

y2=x*exp(x/xw)/(exp(x/xw)+exp(-x/xw));

The quantity xw is the transition width, and has the same units as x. have used this function to model flow through cardiac valves which are basically one-way vales, with ode45(). Choose xw smaller to make the transition sharper.

Plot of y1 and y2 vs. x:

[William Rose]

@Sumit Saha, @Walter Roberson

A smooth symmetric substitute for

sign(x)

is

(exp(x/xw)-exp(-x/xw))/(exp(x/xw)+exp(-x/xw))

Another smooth alternative to sign(x) is

atan(x*pi/(2*xw))*2/pi

which also goes from -1 to +1 and has the same slope at x=0 as the exponential function. xw=transition width (in units of x). Choose xw smaller for a sharper transition. A plot of sign(x) and the smooth functions is below, for xw=1.

The exponential function is better than atan() for this purpose, because the exponential gets to the -1,+1 limits sooner, for the same slope in the center.

[Sumit Saha]

@William Rose Can I use the heavi-side function instead of signum function?

[Walter Roberson]

Can I use the heavi-side function instead of signum function?

You can code whatever you want. Choose the values at random if you want. Choose the values according to the phase of the moon if you want.

However, if your equations are not C2 continuous (continuous second derivatives) then you will get an answer that is incorrect for the situation, because the mathematics behind the ode*() functions assumes C2.

The only time that Heaviside is C2 continuous is if the change between sides is zero.

[William Rose]

@Walter Roberson you make me laugh! We can all use that.

[Walter Roberson]

There turns out to be serious use of choosing values at random, "stochastic differential equations". One of the finance toolboxes (!) has functions for those. But the ode*() routines are not appropriate.

[William Rose]

@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!

Answer 3

[William Rose]

@Sumit Saha

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.

[Sumit Saha]

Thanks a lot...May be ode 23s can work...I'll check that and let you know regarding this issue.