Solving Numerical Integral Implicitly

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

Answer 1

[Walter Roberson]

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.

Answer 2

[yusuf]

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.

[Walter Roberson]

When I substitute in that g, I get the expression evaluating to about -5 * 10^132

It is necessary to use a higher Digits setting as round-off error is crucial in this expression.

The expression has a singularity at g = 225/11 (about 20.54) at which point the denominator goes to 0. My probing appeared to show that the expression as a whole was negative on both sides of 225/11 but looking more closely I see that there is a zero crossing at approximately

20.45391002416497839752127043955070423654750203048625970671478823191317020198663176935360949566221927656712507551824472117488414858009464880684375718019927156153876042420949034484336928538349330556637734407970792045428171493510222258582752843682333646987930

[yusuf]

So, I can use g as a 20.54, right ?

[yusuf]

Also, Walter I need another help from you.

I wrote a code for stationary temperature distribution which is;

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

The T result includes still g , why ? g is 20.54. Because I try to plot T(x). Please help.

Also, Now I try to use perturbation method and find T1 temperature distribution. My equation is;

q*T1''(x)-T1(x)*(f-b*g+i*w*p)=T*(f1-b*g1)-g1

And my prof. say that can be solved by Green's function G(x,y), where the G(x,y) is soluton of;

q*G''(x,y)-G(x,y)*(f-b*g+i*w*p)=Diract(x-y)

And

G(L/2,y)=0 , G(-L/2,y)=0

if I can find the G(x,y), I will get the solution of T1.

Can anyone help me please?

Thank you for everything

Yusuf

[Walter Roberson]

20.45391002416498 will be off by about 1.6/1000 but that's not so bad compared to the slope in that area.

[yusuf]

Can you help me for the other issue ? I wrote above

[Walter Roberson]

It will be later: it is amazing how tiring it is to have a physiotherapist work on your back.

[yusuf]

you should rest good. but if you help me, I will be grateful

[Walter Roberson]

You need to

subs(T,'g',g)

in order to get the T with g replaced.