Solving Numerical Integral Implicitly

7 views (last 30 days)
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

Accepted Answer

Walter Roberson
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.

More Answers (1)

yusuf
yusuf on 4 Feb 2014
The my friend try an another technique , and g is not zero. Can you guys check it ?
f = 0.09;
b = 0.0044;
q = 3.73e-9;
L = 1;
Tw = 250;
Tam = 27;
syms c x g
T = 2*c*cosh(x*((f-b*g)/q)^0.5)+g/(f-b*q);
c = solve(subs(T,'x',L/2)==0,c);
z = simplify(int(subs(T,'c',c),x,-L/2,L/2));
g = solve(z==L*(Tw-Tam),g)
which returns
g =
20.135660961656472105004196502187
We can use double to convert this to floating-point. And I can check that this value of g does indeed solve your equation:
eval(subs(z-L*(Tw-Tam),'g',g)) which returns 0.
  8 Comments
yusuf
yusuf on 4 Feb 2014
you should rest good. but if you help me, I will be grateful
Walter Roberson
Walter Roberson on 10 Feb 2014
You need to
subs(T,'g',g)
in order to get the T with g replaced.

Sign in to comment.

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!