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

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

Answer 1

[Walter Roberson]

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

R1 = 

sR1 = simplify(R1)

sR1 = 

limit(sR1, g, 0)

ans = 

NaN

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

ans = 

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

ans = 

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.

[arvind thakur]

sir how to solve U1

[arvind thakur]

sir i have attached pdf file above can u look at that integration

[Walter Roberson]

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]

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

[arvind thakur]

could any one help me

[Walter Roberson]

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]

sir , due respect i am taking correct but written wrong in code ...i am taking in pascal of stress value and dimension taking in meter .

[arvind thakur]

could anyone plz help me in this

[Walter Roberson]

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]

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]

When I run the calculations through Maple (in case MATLAB is making a mistake), then as soon as I substitute in the specific a and a_0 values, Maple is able to find a closed form integral... and it says that it has a component which is -1i*infinity (assuming that the stress is positive)

[arvind thakur]

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]

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

[arvind thakur]

i am using 2019a

[Walter Roberson]

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]

dear sir i cant find the integral , the paper I have given, how is the result done and the graph is plotted

[arvind thakur]

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

[arvind thakur]

what should i do please give some tips

[Walter Roberson]

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]

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

[Walter Roberson]

u1 = arrayfun(u1fun, X);

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

[arvind thakur]

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

[Walter Roberson]

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:

1. 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
2. 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
3. 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...

1. 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.
2. 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]

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]

You should re-read the paragraph after the P function is defined. It says that for the fibre pressure model, numeric integration was used....sir my question is you have written code by this method only and it will give same result or not? according to that paper plz help sir

[Walter Roberson]

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]

sir due respect plz put the revised and correct code which give the same result so i can take help for further study.

[arvind thakur]

ok it will give same result ?

[Walter Roberson]

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.

[arvind thakur]

what will change in code when I run in matlab 2015 version?

[Walter Roberson]

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]

in 2015 version error came not defined piecewise

[arvind thakur]

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

[arvind thakur]

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]

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

[Walter Roberson]

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

[arvind thakur]

iterative method how to begin if two function dependent on each other , inmy case u(x) and c(x) dependent on each other, so how to integrate u (x)

[Walter Roberson]

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]

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]

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.

Answer 2

[arvind thakur]

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