Solving Numerical Integral Implicitly
7 views (last 30 days)
Show older comments
Hi everyone;
I am working on constant temperature hot-wire anemometry. So I am using second order diff. egn (conduction eqn).
I solved analytically main eqn and found temperature distribution ;
f=0.09;
b=0.0044;
q=3.73E-9;
L=1;
Tw=250;
Tam=27;
T(x)= 2*C1*cosh(x*((f-b*g)/q)^0.5)+g/(f-b*g)
Then C1 has to be determined from boundary condition.
T(+L/2)=0
T(-L/2)=0
Then I found C1 depends on g. Because g is implicitly unknown.
syms c g
solve(2*c*cosh(0.5*(0.09-0.0044*g)/3.73*10^-9)^0.5+g/(0.09-3.73*10^-9*g)==0,c)
g can be determined from constant temperature condition.
1/L*int(T(x)dx,-L/2,L/2)=Tw-Tam
All things considered, my all code is;
clc;
clear all;
f=0.09;
b=0.0044;
q=3.73*10^-9;
L=1;
Tw=250;
Tam=27;
syms c g
c=solve(2*c*cosh(L/2*(0.09-0.0044*g)/3.73*10^-9)^0.5+g/(0.09-0.0044*3.73*10^-9)==0,c)
syms x
z=int(2*c*cosh(x*((f-b*g)/q)^0.5)+g/(f-b*g),x,-L/2,L/2);
g=solve(z==L*(Tw-Tam),g)
g
This condition should give,after performing the integral, an algebraic equation for g. But g (result) is zero . It alyaws gives g as a zero. why ? My Matlab skills are not enough for it. Then I want to plot temperature dist. T(x) . x can be divided 100 parts in length L to plot temperature distribution. Please help for code. Thank you,
Yusuf
0 Comments
Accepted Answer
Walter Roberson
on 3 Feb 2014
The solution for g is not actually 0, but 0 is as close as MATLAB can get to representing it.
If one converts the floating point numbers to rational numbers before using them, then g is approximately
-2.2168581150197961206708392364704289128686142447296*10^(-1062)
The difficulty arises out of your expression
T(x)= 2*C1*cosh(x*((f-b*g)/q)^0.5)+g/(f-b*g)
and in particular the part
cosh(x * sqrt((f-b*g)/q))
when you integrate with x from -L/2 to +L/2 then you end up with the difference between exp(-L*sqrt((f-b*g)/q)) and exp(+L*sqrt((f-b*g)/q)) and that difference needs to balance the other factors of the expression just so. And with those coefficients, the balance is not at 0 but is instead near -2.2E-1062. This is a problem inherent in the expression. e.g.,
int(cosh(x*Y), x = -A .. A)
gives 2*sinh(A*Y)/Y but convert the sinh() to exponential form and you get
(2*((1/2)*exp(A*Y)-(1/2)*exp(-A*Y)))/Y
which is the horrid balancing act I mentioned.
0 Comments
More Answers (1)
See Also
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!