You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Looking to solve this integral with variable bounds
1 view (last 30 days)
Show older comments
%variables
F = .00024; %kg
E = 10^6; %mPa
I = 6.1*10^-9; %kg m^2
L = .01; %m
%solve for Theta sub 0
func = @(OO) L/sqrt((E*I)/(2*F)) - integral(@(theta)1./sqrt(cos(OO)-cos(theta)),OO,pi/2);
g = fzero(func,0);
OS = fsolve(func,.5);
ideally wanted to use fzero but not sure how to do it, the like below, OS = ... kinda works but I keep getting pi/2 and not sure why. Im using this is the follwing equation:
aa = linspace(pi/2*.99, 0, 49);
aa=aa';
%shape of rod given
xi = sqrt((2*E*I)/F)*(sqrt(cos(OS))-sqrt(cos(OS)-cos(aa)));
and not getting the values i am hoping for. any help appreciated
12 Comments
Star Strider
on 11 Jul 2022
I am not certain what you want to do.
In any event, to understand the function, plot it —
F = .00024; %kg
E = 10^6; %mPa
I = 6.1*10^-9; %kg m^2
L = .01; %m
%solve for Theta sub 0
func = @(OO) L/sqrt((E*I)/(2*F)) - integral(@(theta)1./sqrt(cos(OO)-cos(theta)),OO,pi/2);
a = linspace(0, 2*pi, 250);
plotfunc = arrayfun(func, a);
figure
plot(a, real(plotfunc), a, imag(plotfunc))
grid
legend('Real','Imag', 'Location','best')
.
adam puchalski
on 11 Jul 2022
essentially I would like to solve for "OO" in "func" and then use that valuie of OOin my equation for xi.
Torsten
on 11 Jul 2022
Edited: Torsten
on 11 Jul 2022
Up to what you write, you only have to solve for OO once and then insert aa in the equation.
%variables
F = .00024; %kg
E = 10^6; %mPa
I = 6.1*10^-9; %kg m^2
L = .01; %m
%solve for Theta sub 0
func = @(OO) L/sqrt((E*I)/(2*F)) - integral(@(theta)1./sqrt(cos(OO)-cos(theta)),OO,pi/2);
OS = fsolve(func,.5)
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
OS = 1.5708
aa = linspace(0,pi/2*.99, 49);
aa=aa';
%shape of rod given
xi = sqrt((2*E*I)/F)*(sqrt(cos(OS))-sqrt(cos(OS)-cos(aa)));
plot(aa,abs(xi))
Star Strider
on 11 Jul 2022
The two functions you are using for that are root finding algorithms, and the integral function never crosses zero, except at approximately π, as can be seen in the plot.
adam puchalski
on 11 Jul 2022
Edited: Torsten
on 11 Jul 2022
ok so i edited it and ff1 should cross zero now but still does not seem to work
format long
F = .00024; %kg
E = 10^6; %mPa
I = 6.1*10^-9; %kg m^2
L = .01; %m
ep = 10^-8;
z(1) = 0;
for n = 2:10001
z(n) = z(n-1) + ((pi/2)/10000);
ff1(n) = L/sqrt((E*I)/(2*F)) - (2*(ep)^.5)/ (sin(z(n)))^.5 - integral(@(theta) 1./sqrt(cos(z(n))-cos(theta)),z(n)+ep,pi/2);
end
Warning: Minimum step size reached near x = 1.5708. There may be a singularity, or the tolerances may be too tight for this problem.
%figure(5)
plot(z,real(ff1))
func = @(OO) L/sqrt((E*I)/(2*F)) - (2*(ep)^.5)/ (sin(OO))^.5 - integral(@(theta)1./sqrt(cos(OO)-cos(theta)),OO,pi/2);
%OS = fzero(func,.5);
OS = fsolve(func,0.5)
Equation solved, inaccuracy possible.
The vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.
OS =
1.570794630102631
L/sqrt((E*I)/(2*F))
ans =
0.002805147493273
pi/2
ans =
1.570796326794897
adam puchalski
on 11 Jul 2022
the value for OS is either rounded, or inaccurate. Zooming in on the vlaue in the plot, the graph crosses at a value smaller than 1.5708
Torsten
on 11 Jul 2022
Edited: Torsten
on 11 Jul 2022
The graph interpolates because the correct z(n) that makes the function zero is not hit.
I used more points for the graph (see above).
format long
F = .00024; %kg
E = 10^6; %mPa
I = 6.1*10^-9; %kg m^2
L = .01; %m
ep = 10^-8;
func = @(OO) L/sqrt((E*I)/(2*F)) - (2*(ep)^.5)/ (sin(OO))^.5 - integral(@(theta)1./sqrt(cos(OO)-cos(theta)),OO,pi/2);
%OS = fzero(func,.5);
OS = fsolve(func,0.5)
Equation solved, inaccuracy possible.
The vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.
OS =
1.570794630102631
func(OS)
ans =
4.675807340707300e-09
adam puchalski
on 12 Jul 2022
I see, that makes sense! when you fun
func(OS)
thats still not zero. is that due to the step size? can i make the step size even smaller?
Torsten
on 12 Jul 2022
Edited: Torsten
on 12 Jul 2022
It has nothing to do with "step size". At some stage, you reach the limits of floating-point arithmetics. You could try to set TolFun and/or TolX in the options for "fsolve" to a smaller value, but I think that a residual of 5e-9 should suffice for your purposes, doesn't it ?
Be happy that the integrator does not grump because your function has a singularity at "OO" (you try to evaluate 1/sqrt(cos(OO)-cos(OO)) which gives 1/0).
The series expansion around 0, e.g., seems to indicate that your integral does not exist because int(1/x,0,1) = Inf.
syms x
f = 1/sqrt(1-cos(x));
simplify(series(f,x,'ExpansionPoint',0,'order',10))
ans =
Answers (0)
See Also
Tags
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 (한국어)