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
3 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);
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
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 .
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
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
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
dear sir i cant find the integral , the paper I have given, how is the result done and the graph is plotted
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
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
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
sir due respect plz put the revised and correct code which give the same result so i can take help for further study.
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
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
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.
See Also
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 (한국어)