Group Abstract

Message Boards

WOLFRAM COMMUNITY

12.1K Views

3 Replies

0 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Mathematics Mathematica

How to use Nsolve or findroot for a function / over a range

leilazargarzadeh

Posted 12 years ago

Hello Here is the case I don't know how to solve: Over a range of "R", I want to solve etha=Exp[function of etha and R]. That should give an "etha" for each "R" Then calculate some other parameters that are function of both "etha" and "R" I really appreciate any help The code is Clear["Global`"] k = 1.380648810^-23; Av = 6.0210^23; T = 25 + 273.15; Psat1 = 3.16810^3; PL = 101.32510^3; vLinf = (1/997.0479)18.015/1000; gama = 0.07199; theta = Degree 10 ; qc = 1.2210^8; W = 5010^-6; alfa = Degree 45; L = 0.45; Vcontainer = L^3 ; N1 = 3.04910^27; N2 = 3.6710^22; KH = 9.07710^4; bound1 = W/Cos [theta - alfa]; Print["bound1", "=", N[bound1]] (*********************************************************) Vv[R_] := (PiR^3/3)(2 - 3Sin[theta - alfa] + (Sin[ theta - alfa])^3 + (Cos[theta - alfa])^3/Tan[alfa]) Asv[R_] := PiR^2(Cos[theta - alfa])^2/Sin[alfa]; Alv[R_] := 2PiR^2(1 - Sin[theta - alfa]); H[R_] := R(1 - Sin[theta - alfa] + Cos[theta - alfa]/Tan[alfa]); etha[R_] := Exp[(vLinf/(6.0210^23))(PL - Psat1)/(kT) - N2kT/(N1kT - qcVv[R](etha[R]Psat1 - KH))]; N1L[R_] := N1 - qc(etha[R]Psat1Vv[R]/(kT)); N2sL[R_] := (N1L[R]PL/KH); N2L[R_] := N2/(1 + qcPLVv[R]/(N2sL[R]kT)); Rc[R_] := 2gama/(ethaPsat1 + PL(N2L[R]/N2sL[R]) - PL); B[R_] := qc(-2gama * Vv[R]/Rc[R] + gama(Asv[R] Cos [theta] + Alv[R])) + N1kT(N2/N1 - N2L[R]/N1L[R]) + N2kTLog[N1kT/(N1kT - qcethaPsat1Vv[R] + qcVv[R]KH)]; ( I want to plot B vs. H same as Fig 2*) t2 = Table[{ H[R], B[R]/qc}, {R, 10^-6, bound1, 10^-7}]; ListPlot[t2, Joined -> True, PlotStyle -> {Thick}, AxesLabel -> {Style["H", "Graphics", FontSize -> 13], Style["B-B0 /qc (J)", "Graphics", FontSize -> 13]}, AxesOrigin -> {0, 0}, TicksStyle -> Directive[12]]

Hello

Here is the case I don't know how to solve: Over a range of "R", I want to solve etha=Exp[function of etha and R]. That should give an "etha" for each "R" Then calculate some other parameters that are function of both "etha" and "R" I really appreciate any help The code is

Clear["Global`*"]
k = 1.3806488*10^-23;      
Av = 6.02*10^23;                   
T = 25 + 273.15;                       
Psat1 = 3.168*10^3;          
PL = 101.325*10^3;             
vLinf = (1/997.0479)*18.015/1000;   
gama = 0.07199;                      
theta = Degree 10 ;                
qc = 1.22*10^8;
W = 50*10^-6;                        
alfa = Degree 45;
L = 0.45;                                    
Vcontainer = L^3 ;             
N1 = 3.049*10^27;             
N2 = 3.67*10^22;               
KH = 9.077*10^4;        
bound1 = W/Cos [theta - alfa];
Print["bound1", "=", N[bound1]]
(***********************************************************)
Vv[R_] := (Pi*R^3/3)*(2 - 
    3*Sin[theta - alfa] + (Sin[
       theta - alfa])^3 + (Cos[theta - alfa])^3/Tan[alfa])
Asv[R_] := Pi*R^2*(Cos[theta - alfa])^2/Sin[alfa];
Alv[R_] := 2*Pi*R^2*(1 - Sin[theta - alfa]);
H[R_] := R*(1 - Sin[theta - alfa] + Cos[theta - alfa]/Tan[alfa]);

etha[R_] := 
  Exp[(vLinf/(6.02*10^23))*(PL - Psat1)/(k*T) - 
    N2*k*T/(N1*k*T - qc*Vv[R]*(etha[R]*Psat1 - KH))];

N1L[R_] := N1 - qc*(etha[R]*Psat1*Vv[R]/(k*T));     
N2sL[R_] := (N1L[R]*PL/KH);    
N2L[R_] := N2/(1 + qc*PL*Vv[R]/(N2sL[R]*k*T));    

Rc[R_] := 2*gama/(etha*Psat1 + PL*(N2L[R]/N2sL[R]) - PL);    
B[R_] := qc*(-2*gama * Vv[R]/Rc[R] + 
      gama*(Asv[R]* Cos [theta] + Alv[R])) + 
   N1*k*T*(N2/N1 - N2L[R]/N1L[R]) + 
   N2*k*T*Log[N1*k*T/(N1*k*T - qc*etha*Psat1*Vv[R] + qc*Vv[R]*KH)];


(* I want to plot B vs. H same as Fig 2*)
t2 = Table[{ H[R], B[R]/qc}, {R, 10^-6, bound1, 10^-7}];
ListPlot[t2, Joined -> True, PlotStyle -> {Thick}, 
 AxesLabel -> {Style["H", "Graphics", FontSize -> 13], 
   Style["B-B0 /qc (J)", "Graphics", FontSize -> 13]}, 
 AxesOrigin -> {0, 0}, TicksStyle -> Directive[12]]

POSTED BY: leilazargarzadeh

3 Replies

Sort By:

Udo Krause

Udo Krause, NOrganization

Posted 12 years ago

POSTED BY: Udo Krause

leilazargarzadeh

Posted 12 years ago

POSTED BY: leilazargarzadeh

Udo Krause

Udo Krause, NOrganization

Posted 12 years ago

First one has to solve the non-linear functional equation in eta. To avoid the Overflow error message, you define a function In[77]:= Clear[leila] leila[[Eta]_, R_] := Block[{x = (vLinf/Av)(PL - Psat1)/(kT) - N2kT/(N1kT - qcVv[R]([Eta] *Psat1 - KH))}, Exp[x]] this function is rather dull, because you can do In[87]:= Table[{R, FindRoot[[Eta] == leila[[Eta], R], {[Eta], 0}]}, {R, 0, 10, 1}] Out[87]= {{0, {[Eta] -> 1.0007}}, {1, {[Eta] -> 1.00072}}, {2, {[Eta] -> 1.00072}}, {3, {[Eta] -> 1.00072}}, {4, {[Eta] -> 1.00072}}, {5, {[Eta] -> 1.00072}}, {6, {[Eta] -> 1.00072}}, {7, {[Eta] -> 1.00072}}, {8, {[Eta] -> 1.00072}}, {9, {[Eta] -> 1.00072}}, {10, {[Eta] -> 1.00072}}} with In[95]:= leila[1.00072, 8] Out[95]= 1.00072 as well as In[96]:= Table[{R, FindRoot[[Eta] == leila[[Eta], R], {[Eta], 0}]}, {R, 0, 10^8, 10^7}] Out[96]= {{0, {[Eta] -> 1.0007}}, {10000000, {[Eta] -> 1.00072}}, {20000000, {[Eta] -> 1.00072}}, {30000000, {[Eta] -> 1.00072}}, {40000000, {[Eta] -> 1.00072}}, {50000000, {[Eta] -> 1.00072}}, {60000000, {[Eta] -> 1.00072}}, {70000000, {[Eta] -> 1.00072}}, {80000000, {[Eta] -> 1.00072}}, {90000000, {[Eta] -> 1.00072}}, {100000000, {[Eta] -> 1.00072}}} which means Eta is nearly a constant over R with the constants given. Anyway, in case you had a non-constant Eta you could fit its value table {R, Eta} In[97]:= Table[{R, FindRoot[[Eta] == leila[[Eta], R], {[Eta], 0}][[-1, -1]]}, {R, 0,10^8, 10^7}] Out[97]= {{0, 1.0007}, {10000000, 1.00072}, {20000000, 1.00072}, {30000000, 1.00072}, {40000000, 1.00072}, {50000000, 1.00072}, {60000000, 1.00072}, {70000000, 1.00072}, {80000000, 1.00072}, {90000000, 1.00072}, {100000000, 1.00072}} e.g. with Interpolation and use it afterwards because you know it's a function of R now. A symbolic solution for Eta is unavailable as the simplest examples show In[100]:= Reduce[[Eta] == Exp[b/(a - [Eta])], [Eta]] During evaluation of In[100]:= Reduce::nsmet: This system cannot be solved with the methods available to Reduce. >> Out[100]= Reduce[[Eta] == E^(b/(a - [Eta])), [Eta]]

First one has to solve the non-linear functional equation in eta. To avoid the Overflow error message, you define a function

In[77]:= Clear[leila]
leila[[Eta]_, R_] := Block[{x = (vLinf/Av)*(PL - Psat1)/(k*T) - N2*k*T/(N1*k*T - qc*Vv[R]*([Eta] *Psat1 - KH))}, Exp[x]]

this function is rather dull, because you can do

In[87]:= Table[{R, FindRoot[[Eta] == leila[[Eta], R], {[Eta], 0}]}, {R, 0, 10, 1}]
Out[87]= {{0, {[Eta] -> 1.0007}}, {1, {[Eta] -> 1.00072}}, {2, {[Eta] -> 1.00072}}, {3, {[Eta] -> 1.00072}}, 
              {4, {[Eta] -> 1.00072}}, {5, {[Eta] -> 1.00072}}, {6, {[Eta] -> 1.00072}}, {7, {[Eta] -> 1.00072}}, 
              {8, {[Eta] -> 1.00072}}, {9, {[Eta] -> 1.00072}}, {10, {[Eta] -> 1.00072}}}

with

In[95]:= leila[1.00072, 8]
Out[95]= 1.00072

as well as

In[96]:= Table[{R, FindRoot[[Eta] == leila[[Eta], R], {[Eta], 0}]}, {R, 0, 10^8, 10^7}] 
Out[96]= {{0, {[Eta] -> 1.0007}}, {10000000, {[Eta] -> 1.00072}}, {20000000, {[Eta] -> 1.00072}}, 
              {30000000, {[Eta] -> 1.00072}}, {40000000, {[Eta] -> 1.00072}}, {50000000, {[Eta] -> 1.00072}}, 
              {60000000, {[Eta] -> 1.00072}}, {70000000, {[Eta] -> 1.00072}}, {80000000, {[Eta] -> 1.00072}}, 
              {90000000, {[Eta] -> 1.00072}}, {100000000, {[Eta] -> 1.00072}}}

which means Eta is nearly a constant over R with the constants given. Anyway, in case you had a non-constant Eta you could fit its value table {R, Eta}

In[97]:= Table[{R, FindRoot[[Eta] == leila[[Eta], R], {[Eta], 0}][[-1, -1]]}, {R, 0,10^8, 10^7}]
Out[97]= {{0, 1.0007}, {10000000, 1.00072}, {20000000, 1.00072}, {30000000, 1.00072}, 
    {40000000, 1.00072}, {50000000, 1.00072}, {60000000, 1.00072}, {70000000, 1.00072},
     {80000000, 1.00072}, {90000000, 1.00072}, {100000000, 1.00072}}

e.g. with Interpolation and use it afterwards because you know it's a function of R now. A symbolic solution for Eta is unavailable as the simplest examples show

In[100]:= Reduce[[Eta] == Exp[b/(a - [Eta])], [Eta]]
During evaluation of In[100]:= Reduce::nsmet: This system cannot be solved with the methods available to Reduce. >>
Out[100]= Reduce[[Eta] == E^(b/(a - [Eta])), [Eta]]

POSTED BY: Udo Krause

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Feedback