**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# int function return ans in 'int' itself . I got the expression J and L in terms of g but the integration of expression U1 return as in int function. Thanks in advance

8 views (last 30 days)

Show older comments

clc

clear

syms x g

assume(g >= 0)

Ef=427000*10^6; % in MPa

Em=89000*10^6;

Ec=184000*10^6;

mu=0.2825;

ao=0.001;

a=0.00189;

w=0.00512;

x1=0;

%x1=0.001

%R=.0725;

S=198*10^6;

%tow=20;

m1=0.6147+17.1844*(a^2/w^2)+8.7822*(a^6/w^6);

m2=0.2502+3.2889*(a^2/w^2)+70.0444*(a^6/w^6);

%m1=2.7552;

%m2=0.7889;

H=(1+m1*((g-x)/g)+m2*((g-x)/g)^2);

H1=(1+m1*((g-x1)/g)+m2*((g-x1)/g)^2);

% for remote stress S=198 MPa

c=S*((w/(w-ao))+6*w*ao*((0.5*(w-ao)-(x-ao))/(w-ao)^3));

j=int(((S*H)/(sqrt(g-x))),x,[0,ao]);

L=int(((S-c)*H)/(sqrt(g-x)),x, [ao,a]);

R=(H1*j+H1*L);

R1=(R/(sqrt(g-x1)))

U1=int(R1,g,0, a);

##### 0 Comments

### Answers (2)

Walter Roberson
on 25 Sep 2021

Edited: Walter Roberson
on 25 Sep 2021

clc

clear

syms x g

assume(g >= 0)

Ef=427000*10^6; % in MPa

Em=89000*10^6;

Ec=184000*10^6;

mu=0.2825;

ao=0.001;

a=0.00189;

w=0.00512;

x1=0;

%x1=0.001

%R=.0725;

S=198*10^6;

%tow=20;

m1=0.6147+17.1844*(a^2/w^2)+8.7822*(a^6/w^6);

m2=0.2502+3.2889*(a^2/w^2)+70.0444*(a^6/w^6);

%m1=2.7552;

%m2=0.7889;

H=(1+m1*((g-x)/g)+m2*((g-x)/g)^2);

H1=(1+m1*((g-x1)/g)+m2*((g-x1)/g)^2);

% for remote stress S=198 MPa

c=S*((w/(w-ao))+6*w*ao*((0.5*(w-ao)-(x-ao))/(w-ao)^3));

j=int(((S*H)/(sqrt(g-x))),x,[0,ao]);

L=int(((S-c)*H)/(sqrt(g-x)),x, [ao,a]);

R=(H1*j+H1*L);

R1=(R/(sqrt(g-x1)))

sR1 = simplify(R1)

limit(sR1, g, 0)

If you look at the pretty() output, you will see that there is a division by . The numerator has g to a variety of powers. The result is an expression whos limit cannot be taken at 0, and therefore cannot be integrated.

digits(100)

vpa(subs(sR1, g, sym(10)^(-20)))

vpa(subs(sR1, g, sym(10)^(-50)))

You can see from those last couple of evaluations that the real part is steady near 0, but that the imaginary part is exploding. I suspect that the denominator has a root in the range of integration.

##### 40 Comments

Walter Roberson
on 25 Sep 2021

syms x g

assume(g >= 0)

Ef=427000*10^6; % in MPa

Em=89000*10^6;

Ec=184000*10^6;

mu=0.2825;

ao=0.001;

a=0.00189;

w=0.00512;

x1=0;

%x1=0.001

%R=.0725;

S=198*10^6;

%tow=20;

m1=0.6147+17.1844*(a^2/w^2)+8.7822*(a^6/w^6);

m2=0.2502+3.2889*(a^2/w^2)+70.0444*(a^6/w^6);

%m1=2.7552;

%m2=0.7889;

H=(1+m1*((g-x)/g)+m2*((g-x)/g)^2);

H1=(1+m1*((g-x1)/g)+m2*((g-x1)/g)^2);

% for remote stress S=198 MPa

c=S*((w/(w-ao))+6*w*ao*((0.5*(w-ao)-(x-ao))/(w-ao)^3));

j=int(((S*H)/(sqrt(g-x))),x,[0,ao]);

L=int(((S-c)*H)/(sqrt(g-x)),x, [ao,a]);

R=(H1*j+H1*L);

R1=(R/(sqrt(g-x1)))

R1 =

sR1 = simplify(R1)

sR1 =

fplot(real(sR1), [1e-20 1])

title('real part')

G = linspace(1e-20, 1e-2, 250);

R1G = imag(subs(sR1, g, G));

plot(G, R1G, '-*')

title('log imaginary part')

format longg

double(R1G(2:10))

ans = 1×9

1.0e+12 *
2.5261 0.4074 0.1318 0.0556 0.0265 0.0132 0.0064 0.0027 0.0005

double(R1G(end))

ans =

0

idx = find(R1G <= 0, 1)

idx =

11

G(idx), double(R1G(idx))

ans =

0.000401606425702811

ans =

-711872556.23751

vpasolve(imag(sR1), g, 0.0004)

ans =

0.00037621812526461083132509319438441

The imaginary part crosses 0 near 0.000376, but looks to keep going (but might possibly return to 0 at higher values.)

arvind thakur
on 25 Sep 2021

This is the equation I want to integrate for u(x), Have I written the correct code, plz have a look. All the values which are given are in the code.

i want this type of graph and one more thind in u(x) equation ('dl' ) instead of ( 'cl')

I am getting very irritated, plz help me out sir, thanks in advance. in question sigma is S

this is the given value

Walter Roberson
on 27 Sep 2021

Ef=427000*10^6 % in MPa

Ef = 4.2700e+11

If Ef is in megapascals, then

pascals = Ef*1e6

pascals = 4.2700e+17

but according to Table 1,

pascals = 427E9 %427 gigapascals

pascals = 4.2700e+11

You need to decide which units you are using, and you need to make sure that the units are consistent in your calculations.

arvind thakur
on 27 Sep 2021

Walter Roberson
on 29 Sep 2021

The following code takes several minutes to process the release(), as at that point it has try to do the actual integrations. The process is able to get rid of the piecewise(), but it is not at all fast.

My tests so far suggest that the results are guaranteed to have an imaginary part. I am not certain yet if the real part is ever non-zero.

Your version had a lot of confusion between a being a particular value, or a being a symbolic variable. It is necessary for it to be a symbolic variable, as you integrate over L passing in L as the first parameter to H and the first parameter is named a inside of H

syms a a_0 L sigma__inf w x x__prime

assume(a >= 0 & a_0 > 0)

Pi = sym(pi);

a_spm1 = 1.89E-3; % in m %"specimen #1 a = 1.89 mm"

delta_sigma_spm1__inf = 198E6; %in Pa %"specimen #1 delta sigma^infinity = 198 MPa"

a_0_spm1 = 1.0E-3; %in m

w_spm1 = 5.12e-3; % in m

b_spm1 = 2.03E-3; % in m

max_sigma_spm1__inf = 220E6; %in Pa

min_sigma_spm1__inf = 22E6; %in Pa

R_ratio_spm1 = 0.1; % dimensionless

a0_spm2 = 1.0E-3; %in m

w_spm2 = 5.09e-3; % in m

b_spm2 = 1.91E-3; % in m

max_sigma_spm2__inf = 311E6; %in Pa

min_sigma_spm2__inf = 31E6; %in Pa

R_ratio_spm2 = 0.1; % dimensionless

E_f = 427E9; % in Pa

E_m = 89E9; % in Pa

E_c = 184E9; % in Pa

v_c = 0.2825; % dimensionless ratio %was named mu

R = 72.5E-6; % in m %checked

v_f = 0.36; % dimensionless ratio

m1(a) = 0.6147 + 17.1844*(a^2/w^2) + 8.7822*(a^6/w^6); %checked

m2(a) = 0.2502 + 3.2889*(a^2/w^2) + 70.0444*(a^6/w^6); %checked

H(a,x__prime) = 1 + m1(a)*(a-x__prime)/a + m2(a)*(a-x__prime)^2/a^2; %checked

% for remote stress S=198 MPa

c(x) = sigma__inf * 6*w*a_0 * ((w-a_0)/2 - (x - a_0)) / (w - a_0)^3; %checked

P(x__prime) = piecewise(x__prime <= a_0, sigma__inf, sigma__inf - c(x__prime));

u_part2(L) = simplify(int(P(x__prime) * H(L,x__prime) / sqrt(L - x__prime), x__prime, 0, L, 'hold', false))

u_part1(x) = int( H(L,x)/sqrt(L - x) * u_part2(L), L, 0, a, 'hold', true)

u(x) = 2 * (1-v_c^2) / (E_c * Pi) * u_part1(x)

u_spm1(x) = subs(u(x), {a, a_0, sigma__inf, w}, {a_spm1, a_0_spm1, delta_sigma_spm1__inf, w_spm1})

UR1 = release(u_spm1);

%do not generate to file -- R2021b it will generate incorrect code

u1fun = matlabFunction(UR1, 'vars', x);

X = linspace(1e-20, 2e-3, 250);

u1 = arrayfun(u1fun, X);

subplot(2,1,1)

plot(X, real(u1))

title('real part')

subplot(2,1,2)

plot(X, imag(u1))

title('imaginary part')

Walter Roberson
on 29 Sep 2021

Hmmm... Maple says there is a definite integral, and that it involves complex infinity imaginary part.

Your definition for P(x') is troublesome: it does not define a value for x = 0 or for x = a_0 exactly, but the integration ranges look to start from 0.

Walter Roberson
on 29 Sep 2021

arvind thakur
on 30 Sep 2021

your version of code is correct or not sir? it will giving error when i run this code u have written.

Error using sym/int (line 79)

Expected input number 4 to match one of these values:

'IgnoreAnalyticConstraints', 'IgnoreSpecialCases', 'PrincipalValue'

The input, 'hold', did not match any of the valid values.

Error in walter (line 37)

u_part2(L) = simplify(int(P(x__prime) * H(L,x__prime) / sqrt(L - x__prime), x__prime, 0, L, 'hold', false))

Walter Roberson
on 30 Sep 2021

hold and release are documented. Either your MATLAB is corrupted or you are using r2019a or earlier and forgot to mention that. Because of changes to the symbolic toolbox over the years it is important that we know exactly which version you are using.

https://www.mathworks.com/help/symbolic/sym.int.html#mw_d6909eb5-98bd-4510-951b-c9db9323ba60

https://www.mathworks.com/help/symbolic/release.html#mw_d6909eb5-98bd-4510-951b-c9db9323ba60

Walter Roberson
on 30 Sep 2021

syms a a_0 L sigma__inf w x x__prime

assume(a >= 0 & a_0 > 0)

Pi = sym(pi);

a_spm1 = 1.89E-3; % in m %"specimen #1 a = 1.89 mm"

delta_sigma_spm1__inf = 198E6; %in Pa %"specimen #1 delta sigma^infinity = 198 MPa"

a_0_spm1 = 1.0E-3; %in m

w_spm1 = 5.12e-3; % in m

b_spm1 = 2.03E-3; % in m

max_sigma_spm1__inf = 220E6; %in Pa

min_sigma_spm1__inf = 22E6; %in Pa

R_ratio_spm1 = 0.1; % dimensionless

a0_spm2 = 1.0E-3; %in m

w_spm2 = 5.09e-3; % in m

b_spm2 = 1.91E-3; % in m

max_sigma_spm2__inf = 311E6; %in Pa

min_sigma_spm2__inf = 31E6; %in Pa

R_ratio_spm2 = 0.1; % dimensionless

E_f = 427E9; % in Pa

E_m = 89E9; % in Pa

E_c = 184E9; % in Pa

v_c = 0.2825; % dimensionless ratio %was named mu

R = 72.5E-6; % in m %checked

v_f = 0.36; % dimensionless ratio

m1(a) = 0.6147 + 17.1844*(a^2/w^2) + 8.7822*(a^6/w^6); %checked

m2(a) = 0.2502 + 3.2889*(a^2/w^2) + 70.0444*(a^6/w^6); %checked

H(a,x__prime) = 1 + m1(a)*(a-x__prime)/a + m2(a)*(a-x__prime)^2/a^2; %checked

% for remote stress S=198 MPa

c(x) = sigma__inf * 6*w*a_0 * ((w-a_0)/2 - (x - a_0)) / (w - a_0)^3; %checked

P(x__prime) = piecewise(x__prime <= a_0, sigma__inf, sigma__inf - c(x__prime));

u_part2(L) = simplify(int(P(x__prime) * H(L,x__prime) / sqrt(L - x__prime), x__prime, 0, L))

u_part1(x) = int( H(L,x)/sqrt(L - x) * u_part2(L), L, 0, a)

u(x) = 2 * (1-v_c^2) / (E_c * Pi) * u_part1(x)

u_spm1(x) = subs(u(x), {a, a_0, sigma__inf, w}, {a_spm1, a_0_spm1, delta_sigma_spm1__inf, w_spm1})

UR1 = (u_spm1);

%do not generate to file -- R2021b it will generate incorrect code

u1fun = matlabFunction(UR1, 'vars', x);

X = linspace(1e-20, 2e-3, 250);

u1 = arrayfun(u1fun, X);

subplot(2,1,1)

plot(X, real(u1))

title('real part')

subplot(2,1,2)

plot(X, imag(u1))

title('imaginary part')

Not tested in your release

arvind thakur
on 2 Oct 2021

arvind thakur
on 2 Oct 2021

also i have given screenshot and i am attaching paper in which that integration is available

Walter Roberson
on 2 Oct 2021

You should re-read the paragraph after the P function is defined. It says that for the fibre pressure model, numeric integration was used. That means that they do not have a closed form integral.

The paragraph goes on to say that for the other model they had to either use an iterative approach or else a finite element approach. Both approaches gave similar answers, but finite element was faster so they used that. Again this means that they do not have a closed form solution.

arvind thakur
on 2 Oct 2021

how to do numerical integration for this i am new on matlab plz help me sir

Walter Roberson
on 2 Oct 2021

u1 = arrayfun(u1fun, X);

The result of that call is the numeric integration evaluated at all the points in X.

arvind thakur
on 2 Oct 2021

It's taking to much time on my laptop sir so any modification so time duration will reduce.?

Walter Roberson
on 2 Oct 2021

Unfortunately, MATLAB is very slow to examine the outer symbolic integral and decide that it does not know how to solve it.

MATLAB has a few tricks available to postpone evaluation of a symbolic integral:

- Since R2016b, vpaintegral() can be used. vpaintegral() will examine the expression, and if it sees any unbound symbolic variables other than the variable of integration, then it will leave the integration unevaluated. In the case of having two levels of integration, if you have three total symbolic variables (two being integrated over plus one other) then it would not try to calculate the outer integral -- not until you substituted in a specific number value for the third variable. When the integration does eventually go ahead, it will be done using symbolic floating point numbers
- Since R2019b, when you use int() you can specify 'hold', true and MATLAB will postpone integration until you use release() to ask that it go ahead with the integration. When the integration does eventually go ahead, it will be done using full indefinite precision, looking for a closed form formula
- if symbolic integration using int() decides that the integral is too complicated for it, then it will return the unevaluated integral in the form of a symbolic int(). With some struggle, you can force a complicated-but-feasible symbolic integral to not be feasible so that symbolic int() is returned (probably after a long time), and then later force that term to drop out... but it probably will then waste time doing the integral all over again, so you mostly do not gain much. And the calculation would be done in full symbolic indefinitely precise mode

Ideally for performance you would like to be able to take the symbolic integral form and use matlabFunction() to convert it to an anonymous function that does numeric integration. However...

- matlabFunction() looks at the vpaintegral() call and gives an error saying it does not know what the "_range" operation is. "_range" is the internal way vpaintegral records the limit of integration. I don't know why Mathworks did not make this work years ago.
- matlabFunction() looks at the 'hold' option of int() and says it doesn't know how to do that

So all you are left with, in full symbolic form, is to allow the int() to go ahead without vpaintegral() and without 'hold', and waste a lot of time (possibly hours) figuring out that it doesn't know how to do the integration... after which matlabFunction() is then willing to convert the int() into a call to integral() . vectorization is not handled, so you have to take care to arrayfun() the call to the anonymous function, but you can get something . But only after you have wasted a lot of time while int() figures out that no, the integral really is too complicated for it.

Walter Roberson
on 2 Oct 2021

So what are your options for performance?

Well, it turns out that the inner integral u_part2 can be done without much trouble. So you ask for that to be done. And then you matlabFunction() that, and hand wrap it inside an anonymous function that uses integral()

Or...

You write the whole thing in terms of anonymous functions and integral() calling integral() to start with, skipping any symbolic step.

At the moment I have my system calculating the

u(x) = 2 * (1-v_c^2) / (E_c * Pi) * u_part1(x)

step. When it finishes, I will attach the result as a .mat file.

I will also have it do

u_spm1(x) = subs(u(x), {a, a_0, sigma__inf, w}, {a_spm1, a_0_spm1, delta_sigma_spm1__inf, w_spm1})

to substitute in the parameters for "speciman 1" and I will attach that too.

Once you have u_spm1() (which does the symbolic calculation using the parameters for speciman 1), evaluation using matlabFunction() is not exactly "fast", but it is much much faster than going the vpaintegral() route.

arvind thakur
on 4 Oct 2021

Walter Roberson
on 4 Oct 2021

The code I already posted does numeric integration. https://www.mathworks.com/matlabcentral/answers/1460464-int-function-return-ans-in-int-itself-i-got-the-expression-j-and-l-in-terms-of-g-but-the-integra#comment_1761304

The matlabFunction() call converts the symbolic integration calls into numeric integration calls.

arvind thakur
on 4 Oct 2021

Walter Roberson
on 4 Oct 2021

You plan to use floating point numbers. When you use floating point numbers,

a + b + c

does not give the same result as

a + c + b

so you should not expect that two different floating point implementations will give the same result.

Walter Roberson
on 4 Oct 2021

Nothing. I already posted the code that is compatible with R2015a.

The version that used "hold" needed R2019b or later, but I already took that into account when I posted https://www.mathworks.com/matlabcentral/answers/1460464-int-function-return-ans-in-int-itself-i-got-the-expression-j-and-l-in-terms-of-g-but-the-integra#comment_1761304

R2015a might not simplify the code as well, so it might produce a different result due to floating point round-off.

arvind thakur
on 6 Oct 2021

any update sir, in my laptop the code is not giving result its taking too much time.

arvind thakur
on 6 Oct 2021

sir if I am doing integration as indefinte integral its give answer in the form of variable

so how we put value of variable at different stages.

arvind thakur
on 8 Oct 2021

sir how to do this coding with shear lag model could plz suggest to begin

Walter Roberson
on 8 Oct 2021

Edited: Walter Roberson
on 9 Oct 2021

I have never looked at shear log models, so I am not able to make any suggestions for that.

arvind thakur
on 8 Oct 2021

Walter Roberson
on 9 Oct 2021

I have done some more testing...

In R2021b, the code version I posted in https://www.mathworks.com/matlabcentral/answers/1460464-int-function-return-ans-in-int-itself-i-got-the-expression-j-and-l-in-terms-of-g-but-the-integra#comment_1760214 that uses the 'hold' option, then substitutes in numeric values, and then does release(), is able to do the integral in about 2/3 of an hour.

However, the code version I posted in https://www.mathworks.com/matlabcentral/answers/1460464-int-function-return-ans-in-int-itself-i-got-the-expression-j-and-l-in-terms-of-g-but-the-integra#comment_1761304 which is the same except not using the 'hold' option (because your version does not support it), takes... days. Long enough that I had to kill the processing to install an operating system security update.

I do not think you are going to be able to use that approach in your older version. Possibly if you substitute in the numeric values before trying to do the integrations it might work for you in a couple of hours.

Walter Roberson
on 9 Oct 2021

Before R2016b, to do piecewise() you can use

symLIST = @(varargin)feval(symengine,'DOM_LIST',varargin{:});

and then

feval(symengine, 'piecewise', symLIST(Condition1,Result1), symLIST(Condition2,Result2))

Walter Roberson
on 9 Oct 2021

inmy case u(x) and c(x) dependent on each other

No they do not. u(x) depends upon c(x) but c(x) does not depend upon u(x) . An iterative approach is not called for.

arvind thakur
on 9 Oct 2021

here is the c(x), ignore previous c(x) and all other value and equation is same.

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)