Problem with code - symbolic integration takes awfully long time

Hi all, I am facing a problem concerning symbolic integration.

I have a function that I need to double integrate for variables x1 and y1. This is what my function looks like:

fun_o = p_outer*mi*(1-exp(-jo/K))*sin_delta1;
ymin = -L/2+cy-so;
ymax = L/2+cy-so

;

unknown (symbolic) parameters are x1, y1 and so and:

p_outer - is a function of y1 and so

jo - is a function of x1, y1 and so

sin_delta1 - is a function of y1

so in not dependent on neither x1 nor y1. Other parameters are numerical values.

This is what my integral looks like:

Fy_outer = int(int(fun_o, x1, -b/2, b/2), y1, ymin, ymax);

After I run the code, it takes awfully long time to compute this integral. Plus I need three more integrals which are as complex as this or even more complex than this integral. Plus I would need to run the code multiple times with different numerical values.

Profiler shows that certain mupadex function consumes the time (the 400s is the most I have let it run):

Is there any way to speed up the simulation, what am I doing wrong?

[Walter Roberson]

To confirm, you have this structure

syms x1 y1 so real
syms b positive
syms cy K L mi real
syms p_outer(y1, so)
syms jo(x1, y1, so)
syms sin_delta1(y1)
fun_o = p_outer(y1, so)*mi*(1-exp(-jo(x1,y1,so)/K))*sin_delta1(y1)

fun_o = 

ymin = -L/2+cy-so

ymin = 

ymax = L/2+cy-so

ymax = 

Fy_outer = int(int(fun_o, x1, -b/2, b/2), y1, ymin, ymax)

Fy_outer = 

except with known sin_delta, p_outer, and jo ?

We are going to need to know the functions to test further; also knowing the restrictions on the variables such as cy and K would help (for example can we assume that L > 0 ?)

[Ste_M]

My 150-line code was working perfectly fine, and executing everything in about 5 seconds until I added this double integral so I did not want to burden the question with additional code parts that work.

But for better understanding I will add everything related to this integral (be aware that parameter names are a little bulky).

Input: (note that x2 and y2 are related to the second part of code which is yet to be written);

syms so Omega_outer x1 y1 t Omega_inner x2 y2
mi = 0.9;
K = 0.075;
B = 2.54;
L = 3.8;
V = 1; 
R = 15; 
W = 25500;
b = 0.45;
cy = 0;
cx = 0;
h = 1.3;
fr = 0.0263;
r = 0.32;
g = 9.981;

Function is:

fun_o = p_outer*mi*(1-exp(-jo/K))*sin_delta1;

Brakedown of the function:

• Pressure distribution function is as follows:

p_outer = W_outer / (b*L) - (12/(b*L^3))*(W/2*cy+h*W*V^2*sin(Beta)/(2*g*R))*(y1+so-cy);

So p_outer is a function of so and y1 because Beta=f(so), or to be precise:

Beta = asin(so/R);

• mi is a number
• jo is a function of x1, so and y1 (note that Omega_outer is not dependent on either variable - acts as constant, but is a symbolic parameter) - OmegaR = numerical value = V/R

R_prim = R*cos(Beta);
jXo = (R_prim + B/2 + cx + x1) * (cos((L/2 + cx - so - y1)*OmegaR/(r*Omega_outer))-1) - y1*sin((L/2 + cy - so - y1)*OmegaR/(r*Omega_outer));
jYo = (R_prim + B/2 + cx + x1) * (sin((L/2 + cy - so - y1)*OmegaR/(r*Omega_outer))) - (L/2 + cy - so) + y1*cos((L/2 + cy - so - y1)*OmegaR/(r*Omega_outer));

jo = sqrt(jXo^2 + jYo^2);

• sin_delta1 is function of x1, y1 and so (my apologies, I made a mistake in original question, I forgot R_prim which is function of so)

Vjyo = (R_prim + B/2 + cx + x1)* OmegaR - r*Omega_outer;
Vjxo = -y1*OmegaR;

sin_delta1 = Vjyo/sqrt(Vjxo^2 + Vjyo^2);

• K is a numerical value

To sum up - symbolic parameters are x1, y1, so, Omega_outer

- everything else is a numerical value

When I stop execution it always shows this

Operation terminated by user during sym/int (line 162)

In Wong_2ndtry (line 161)

Fy_outer = int(int(fun_o, x1, -b/2, b/2), y1, ymin, ymax);

So I am pretty sure that the integral is the problem.

p_outer could be simplified to be uniformly distributed, as opposed to current trapezoidal form. Integral2 can not be used as it requires numerical values. If needed I am willing to provide full code.

Thanks in advance!

[Walter Roberson]

Sometimes you can improve performance by asking to delay evaluation of the inner integral

syms so Omega_outer x1 y1 t Omega_inner x2 y2
mi = 0.9;
K = 0.075;
B = 2.54;
L = 3.8;
V = 1; 
R = 15; 
W = 25500;
b = 0.45;
cy = 0;
cx = 0;
h = 1.3;
fr = 0.0263;
r = 0.32;
g = 9.981;

syms OmegaR W_outer real

Beta = asin(so/R);

p_outer = W_outer / (b*L) - (12/(b*L^3))*(W/2*cy+h*W*V^2*sin(Beta)/(2*g*R))*(y1+so-cy);

R_prim = R*cos(Beta);
jXo = (R_prim + B/2 + cx + x1) * (cos((L/2 + cx - so - y1)*OmegaR/(r*Omega_outer))-1) - y1*sin((L/2 + cy - so - y1)*OmegaR/(r*Omega_outer));
jYo = (R_prim + B/2 + cx + x1) * (sin((L/2 + cy - so - y1)*OmegaR/(r*Omega_outer))) - (L/2 + cy - so) + y1*cos((L/2 + cy - so - y1)*OmegaR/(r*Omega_outer));

jo = sqrt(jXo^2 + jYo^2);

Vjyo = (R_prim + B/2 + cx + x1)* OmegaR - r*Omega_outer;
Vjxo = -y1*OmegaR;

sin_delta1 = Vjyo/sqrt(Vjxo^2 + Vjyo^2);

fun_o = p_outer*mi*(1-exp(-jo/K))*sin_delta1

fun_o = 

Fy_inner = int(fun_o, x1, -b/2, b/2, 'hold', true)

Fy_inner = 

ymin = -L/2+cy-so

ymin = 

ymax = L/2+cy-so

ymax = 

Fy_outer = int(Fy_inner, y1, ymin, ymax, 'hold', true)

Fy_outer = 

result = release(Fy_outer)

I do not show the result of the release here because it would exceed the time limit.

[Walter Roberson]

On my system, the release took about 216 seconds, and did not change the value. There is a chance that if OmegaR and W_outer were defined with numeric values that it could get further.

... though when I look at the denominator, it seems unlikely to me that you would be able to get anywhere symbolically even if numeric values for those were known.

It might also be worth experimenting with reversing the order of integration.

It looks to me as if your best chance at a symbolic solution would be if so were defined numerically as well as OmegaR

[Ste_M]

OmegaR is defined numerically as V/R which is 1/15. Defining so is the ultimate solution of the problem. I am trying to formulate a code to evaluate contact forces between two bodies. The three variables unknown being so, Omega_outer and Omega_inner which can be obtained from the equilibrium equations of motion. Problem is that I need Fy_outer and three more similar formulations (Fy_inner and two moments which are basically double integrals of Fy_outer*x1 and Fy_inner*x2) in order to set up system of three equations. All of these formulations depend on so, so I need to express them in the form f=f(so), put them in the system of equations, define so and then evaluate forces when so is known for the given set of inputs.