I can't integrate symbolic expression with abs(x), and using sign(real(x)) produces wrong results

I am using symbolic expression and I need to integrate a function that contains abs(), this takes so long without producing any result, I tried using sign(x) instead I have same result as using abs(), I tried using sign(real(x)) and this produces results but with error in the results.

Please, can anyone help with this?

Result when I use sign(real(x))

Answer 1

[Star Strider]

My version of MATLAB cannot run images of code.

Apparently the variables have numeric values or the plot would not exist, so one option is to use matlabFunction on ‘drag_cw_function_column’ and integrate the resulting anonymous function numerically.

.

[Sulaiman Ogungboyega]

@Star Strider, Yes the variables have a numeric values, but this function is just a part of other functions that adds up to the final plot, and I have to plot how the final_function (containing drag and other functions) changes over the variables.

So, I will need to keep them as variables till I can subs the values in the final_function.

[Walter Roberson]

Please post the code as text so we can test.

[Star Strider]

There are other problems, most notably that ‘drag_cw_function_column’ is a function fo three variables, being integrated with respect only to ‘z’.

What are the plans for the rest of the variables, and what is the desired end result?

[Sulaiman Ogungboyega]

@Star Strider the plan is subs the value for the variables when I have the final function. I need to check how the variables affect the result for the final_function, not just how they affect drag.

I am attaching a part of the code as "help.txt" the same trend can be seen on the fplot.

This is just a part of the code, I really need to keep them as variables till I get to final function or I will really have a lot of loops which will make the code very cumbersome.

Note: I noticed if the values are changed (amplitude set to 3 for instance) the plots come out good.

[Star Strider]

I would appreciate a bit more information about this, since I can’t understand what the code does.

I looked at several sub-expressions in the function with the discontinuities, and have so far been unable to find the source.

rho = 1030; 
C_d = 0.7;
C_m = 1.0; 
omega = pi/6;
k_cw= 0.0242;
U_current = 1.5;
amplitude = 7.5;
diameter = [6, 8];
depth = 350;
phase_velocity_cw = 21.6184;
wave_length_cw = 259.4214;

du_dt_cw_max = @(z) amplitude*omega^2*(cosh(k_cw*(depth+z)))*(1-(U_current/phase_velocity_cw))/sinh(k_cw*depth);
inertia_cw_function_column = @(z) 0.25*C_m*rho*pi*diameter(1).^2*du_dt_cw_max(z);

syms z x t

% Forces on the column 
phase_cw(x,t) =  omega*(t) -k_cw*(x);

u_wave_current (z,x,t) = U_current+ (amplitude*omega*(cosh(k_cw*(depth+z)))*(1-(U_current/phase_velocity_cw))/sinh(k_cw*depth)).*(sin(phase_cw(x,t)));
drag_cw_function_column_trial(z,x,t)= 0.5*C_d*rho*diameter(1)*(u_wave_current(z,x,t)).^2*sign(real(u_wave_current(z,x,t)));
drag_cw_function_column = matlabFunction(drag_cw_function_column_trial);

figure ("Name","Plot of drag to chech sign(u)")
fsurf(drag_cw_function_column(z,x,0),[0 200 -34 0]); % plot of drag to show sign(u) is effective and drag is negative where needed
title ("Plot of drag to chech sign(u)");

%%
F_drag_column_cw = @(z,x,t) drag_cw_function_column(z,x,t); %general expression for F_drag on column at any point per unit length

F_inertia_column_cw = @(z,x,t) inertia_cw_function_column (z)*cos(phase_cw(x,t)); %general expression for F_inertial on column at any point
F_horizontal_cw_column = @(z,x,t) F_drag_column_cw(z,x,t) + F_inertia_column_cw(z,x,t); %general expression for Total horizontal force at any point on the column at any point per unit length 

F_horizontal_moment(z,x,t) = F_horizontal_cw_column (z,x,t)*(z+34);
int_moment_horizontal_column(z,x,t) = int(F_horizontal_moment,z); %intergrated general expression over z
sum_horizontal_moment_column(x,t) = subs(int_moment_horizontal_column,z,0)-subs(int_moment_horizontal_column,z,-30); %integrated over interval -30 to 0

Moment_horizontal_column(x) = sum_horizontal_moment_column(x,0);

figure("Name","plot showing similar trend")
fplot(Moment_horizontal_column(x),[0 wave_length_cw]); % the plot shows similar trend like the surf plot posted earlier 
title("Plot showing similar behaviour")

.

[Sulaiman Ogungboyega]

@Star Strider, the code caluclates the moment on a column across the wave length of a wave.

Moment is a multiple of effective distance and the force on the column, which is a sum of drag force and inertia force on the column.

I hope this helps.

[Sulaiman Ogungboyega]

The source of the discontinuity is needing to integrate expression "drag_cw_function_column_trial" containing sign(real(u_wave_current(z,x,t))). If that part of the expression was removed, the code runs fine without discontinuity.

[Star Strider]

I still don’t understand the code or the arguments.

Simply using ‘real(u_wave_current(z,x,t))’ produces a smooth sinusiodal curve with no discontinuities that the sign call introduces. Is there anything wrong with just using that instead of thresholding it? Also, that function appears to me to be pure real (although I have no idea what the values of ‘z’ and ‘t’ are supposed to be, so it could be complex in some instances).

rho = 1030; 
C_d = 0.7;
C_m = 1.0; 
omega = pi/6;
k_cw= 0.0242;
U_current = 1.5;
amplitude = 7.5;
diameter = [6, 8];
depth = 350;
phase_velocity_cw = 21.6184;
wave_length_cw = 259.4214;

du_dt_cw_max = @(z) amplitude*omega^2*(cosh(k_cw*(depth+z)))*(1-(U_current/phase_velocity_cw))/sinh(k_cw*depth);
inertia_cw_function_column = @(z) 0.25*C_m*rho*pi*diameter(1).^2*du_dt_cw_max(z);

syms z x t

% Forces on the column 
phase_cw(x,t) =  omega*(t) -k_cw*(x);

u_wave_current (z,x,t) = U_current+ (amplitude*omega*(cosh(k_cw*(depth+z)))*(1-(U_current/phase_velocity_cw))/sinh(k_cw*depth)).*(sin(phase_cw(x,t)));
drag_cw_function_column_trial(z,x,t)= 0.5*C_d*rho*diameter(1)*(u_wave_current(z,x,t)).^2*sign(real(u_wave_current(z,x,t)));
drag_cw_function_column = matlabFunction(drag_cw_function_column_trial);

figure ("Name","Plot of drag to chech sign(u)")
fsurf(drag_cw_function_column(z,x,0),[0 200 -34 0]); % plot of drag to show sign(u) is effective and drag is negative where needed
title ("Plot of drag to chech sign(u)");

%%
F_drag_column_cw = @(z,x,t) drag_cw_function_column(z,x,t); %general expression for F_drag on column at any point per unit length

F_inertia_column_cw = @(z,x,t) inertia_cw_function_column (z)*cos(phase_cw(x,t)); %general expression for F_inertial on column at any point
F_horizontal_cw_column = @(z,x,t) F_drag_column_cw(z,x,t) + F_inertia_column_cw(z,x,t); %general expression for Total horizontal force at any point on the column at any point per unit length 

F_horizontal_moment(z,x,t) = F_horizontal_cw_column (z,x,t)*(z+34);
int_moment_horizontal_column(z,x,t) = int(F_horizontal_moment,z); %intergrated general expression over z
sum_horizontal_moment_column(x,t) = subs(int_moment_horizontal_column,z,0)-subs(int_moment_horizontal_column,z,-30); %integrated over interval -30 to 0

Moment_horizontal_column(x) = sum_horizontal_moment_column(x,0);

figure("Name","plot showing similar trend")
fplot(Moment_horizontal_column(x),[0 wave_length_cw]); % the plot shows similar trend like the surf plot posted earlier 
title("Plot showing similar behaviour")

z = 1;
t = 1;
figure
fplot(sign(real(u_wave_current(z,x,t))), [0 wave_length_cw])
grid
ylim([-1 1]*1.1)

figure
fplot((real(u_wave_current(z,x,t))), [0 wave_length_cw], '-b')
hold on
% fplot((imag(u_wave_current(z,x,t))), [0 wave_length_cw], '-g')
fplot((abs(u_wave_current(z,x,t))), [0 wave_length_cw], '-r')
hold off
grid

.

[Sulaiman Ogungboyega]

So basically I need an expression that squares u(z,x,t) and retains the sign, such that if U(z,x,t) = -2, the expression performs u^2 and returns -4, not 4. If U(z,x,t) is 3 u^2 returns 9. I need to also be able to perform int() on the expression.

Using u*abs(u) didn’t work as when I try to int() the expression containing u*abs(), Matlab does not run it.

Sign(U(z,x,t)) only returns a +1 or -1 depending on the value on U, so I am using u^2*sign(u) in place of u*abs(u).

I tried using sign(real(u)) because Matlab took so long running just sign(u), the value of U will always be real.

I hope this helps.

[Walter Roberson]

syms v
u = cos(v) - sin(v)   %just some mathematical expression

u = 

y = u*abs(u)

y = 

int(y, v, [-3 3])

ans = 

MATLAB is able to integrate u*abs(u) when I try, for simple variables u, or for a non-trivial but not overly-complex u ?

[Sulaiman Ogungboyega]

I am not sure, I used u*abs(u) in the "drag_cw_function_trial" and matlab just keeps running without any result. You should be able to reproduce same thing in the code by changing the "drag_cw_function_trial"

[Walter Roberson]

td = tempdir();
trialfile = fullfile(td, 'drag_cw_function_column_trial.m');
addpath(td);

rho = 1030; 
C_d = 0.7;
C_m = 1.0; 
omega = pi/6;
k_cw= 0.0242;
U_current = 1.5;
amplitude = 7.5;
diameter = [6, 8];
depth = 350;
phase_velocity_cw = 21.6184;
wave_length_cw = 259.4214;

du_dt_cw_max = @(z) amplitude*omega^2*(cosh(k_cw*(depth+z)))*(1-(U_current/phase_velocity_cw))/sinh(k_cw*depth);
inertia_cw_function_column = @(z) 0.25*C_m*rho*pi*diameter(1).^2*du_dt_cw_max(z);

syms z x t

% Forces on the column 
phase_cw(x,t) =  omega*(t) -k_cw*(x);

u_wave_current (z,x,t) = U_current+ (amplitude*omega*(cosh(k_cw*(depth+z)))*(1-(U_current/phase_velocity_cw))/sinh(k_cw*depth)).*(sin(phase_cw(x,t)));
drag_cw_function_column_trial(z,x,t)= 0.5*C_d*rho*diameter(1)*(u_wave_current(z,x,t)).^2*piecewise(real(u_wave_current(z,x,t))<0,-1,1);
drag_cw_function_column = matlabFunction(drag_cw_function_column_trial, 'File', trialfile, 'Optimize', false);

figure ("Name","Plot of drag to chech sign(u)")
fsurf(drag_cw_function_column_trial(z,x,0),[0 200 -34 0]); % plot of drag to show sign(u) is effective and drag is negative where needed
title ("Plot of drag to chech sign(u)");

%%
F_drag_column_cw = drag_cw_function_column; %general expression for F_drag on column at any point per unit length

F_inertia_column_cw = @(z,x,t) inertia_cw_function_column (z)*cos(phase_cw(x,t)); %general expression for F_inertial on column at any point
F_horizontal_cw_column = @(z,x,t) drag_cw_function_column_trial(z,x,t) + F_inertia_column_cw(z,x,t); %general expression for Total horizontal force at any point on the column at any point per unit length 

F_horizontal_moment(z,x,t) = F_horizontal_cw_column (z,x,t)*(z+34);
int_moment_horizontal_column(z,x,t) = int(F_horizontal_moment,z) %intergrated general expression over z

int_moment_horizontal_column(z, x, t) = 

sum_horizontal_moment_column(x,t) = subs(int_moment_horizontal_column,z,0)-subs(int_moment_horizontal_column,z,-30) %integrated over interval -30 to 0

sum_horizontal_moment_column(x, t) = 

Moment_horizontal_column(x) = sum_horizontal_moment_column(x,0)

Moment_horizontal_column(x) = 

figure("Name","plot showing similar trend")
fplot(Moment_horizontal_column(x),[0 wave_length_cw]); % the plot shows similar trend like the surf plot posted earlier 
title("Plot showing similar behaviour")

z = 1;
t = 1;
figure
fplot(sign(real(u_wave_current(z,x,t))), [0 wave_length_cw])
grid
ylim([-1 1]*1.1)

figure
fplot((real(u_wave_current(z,x,t))), [0 wave_length_cw], '-b')
hold on
% fplot((imag(u_wave_current(z,x,t))), [0 wave_length_cw], '-g')
fplot((abs(u_wave_current(z,x,t))), [0 wave_length_cw], '-r')
hold off
grid

[Walter Roberson]

You can see from the output then when I used piecewise(), that all of the int() were able to succeed. The result is piecewise, but nothing overly complex.

[Star Strider]

@Walter Roberson — I’ve been a bit lost here for the last day or so because I don’t understand what the code does, only that thus far I’ve not been able to make any headway (significantly, at least) in solving it. So I appreciate your insights.

[Walter Roberson]

@Star Strider

I do not know what the code does either ;-) I am just here for the part about abs() potentially causing problems in integration.

[Star Strider]

@Walter Roberson — I struggled with that as well, especially since the part of the code that seems to be causing the problem appears to be purely real so the abs() simply creates a discontinuity. Adding 2 to it (then perhaps scaling it if the amplitude is critical so that it is scaled to span (0,5)) would make it positive with no discontinuities. I cannot figure out what’s causing the other discontinuities. I at first thought that might be some sort of ‘wrapping¹ effect of a trigonometric function, however that appears to be unlikely.

[Sulaiman Ogungboyega]

Thank you @Star Strider and @Walter Roberson, I had to go with numerical expression for this.

[Star Strider]

I figured that was the best option.

(I still don’t understand the code!)