Thanks Preben

Amir

I get the correct results.

Thanks for your time and efforts for my code.

Amir

Dear Proben

I get this from you response. but this question is different a little from that I ask you in contact. Here, I want to modify the code in a way that phi0 is calculated with the following addition constraint

evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet))) / evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet)))-0.02=0

Also, I use the secant method becuase it is very quick and fsolve is time consuming.

Thanks

Amir

 Wow!

It was very interesting. It is very quick.

Thanks a lot.

Amir

 Wow!

It was very interesting. It is very quick.

Thanks a lot.

Amir

@Carl Love 

I'm agree with you

@Carl Love 

I'm agree with you

Hi

thanks for your previous comments

however, I want to add a more restriction condition to the above code. In this way, I want to find the value of phi(0) in such a way that evalf(2/(1-zet^2)*(Int((1-eta)*F1(eta),eta=0..1-zet)))-0.02=0. how can i add this?

thanks for your attention in advance

Hi

thanks for your previous comments

however, I want to add a more restriction condition to the above code. In this way, I want to find the value of phi(0) in such a way that evalf(2/(1-zet^2)*(Int((1-eta)*F1(eta),eta=0..1-zet)))-0.02=0. how can i add this?

thanks for your attention in advance

@Carl Love 

yes. it is the same problem. actually the current procedure works. thanks to you and preben.

@Carl Love 

yes. it is the same problem. actually the current procedure works. thanks to you and preben.

Dear Preben

Thanks. it works

Amir

Dear Preben

Thanks. it works

Amir

Dear Carl

thanks for your reply

actually I can not understand your code appropraitely. (it is advanced)

But I Wrote a Code and found a solution to solve these equations.

But it doesnt find the _J_J_ automatically and I correct it manually.

the following code is my originla code

restart:
Digits := 12;
nta:=20;
gama1:=0:
phi0:=0.04;
                               12
                               20
                              0.04

rhocu:=2/(1-zet^2)*int((1-eta)*rho(eta)*c(eta)*u(eta),eta=0..1-zet);
#ff:=unapply((1-eta)*rho(eta)*c(eta)*u(eta),eta);
#rhocu:=2/(1-zet^2)*(1-zet)*( ff(0)/2 +sum(ff(0+ik*(1-zet)/nta),ik=1..nta-1))/nta;


eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])+((1/(eta-1)+1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
eq2:=diff(T(eta),eta,eta)+1/(k(eta)/k1[w])*(2/(1-zet^2)*rho(eta)*c(eta)*u(eta)/p2+( (a[k1]+2*b[k1]*phi(eta))/(1+a[k1]*phi1[w]+b[k1]*phi1[w]^2)*diff(phi(eta),eta)-k(eta)/k1[w]/(1-eta)*diff(T(eta),eta) ));
eq3:=diff(phi(eta),eta)-phi(eta)/(N[bt]*(1-gama1*T(eta))^2)*diff(T(eta),eta);
 2 (int((1 - eta) rho(eta) c(eta) u(eta), eta = 0 .. 1 - zet))
 -------------------------------------------------------------
                                  2                           
                           1 - zet                            
    /  d   /  d         \\   mu1[w]
    |----- |----- u(eta)|| + -------
    \ deta \ deta       //   mu(eta)

         /                 /  d           \\               
         |          mu_phi |----- phi(eta)||               
         |   1             \ deta         /| /  d         \
       + |------- + -----------------------| |----- u(eta)|
         \eta - 1           mu(eta)        / \ deta       /
                                /      /                        
                                |      |                        
/  d   /  d         \\     1    |      |2 rho(eta) c(eta) u(eta)
|----- |----- T(eta)|| + ------ |k1[w] |------------------------
\ deta \ deta       //   k(eta) |      |     /       2\         
                                \      \     \1 - zet / p2      

                                /  d           \
     (a[k1] + 2 b[k1] phi(eta)) |----- phi(eta)|
                                \ deta         /
   + -------------------------------------------
                                          2     
         1 + a[k1] phi1[w] + b[k1] phi1[w]      

            /  d         \\\
     k(eta) |----- T(eta)|||
            \ deta       /||
   - ---------------------||
        k1[w] (1 - eta)   ||
                          //
                                       /  d         \
                              phi(eta) |----- T(eta)|
           /  d           \            \ deta       /
           |----- phi(eta)| - -----------------------
           \ deta         /            N[bt]         
mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta);
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta);
rhop:=3880;
rhobf:=998.2;
cp:=773;
cbf:=4182;
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta);
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta);
mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta));

a[mu1]:=39.11;
b[mu1]:=533.9;
mu1[bf]:=9.93/10000;
a[k1]:=7.47;
b[k1]:=0;
k1[bf]:=0.597;
zet:=0.5;
#phi(0):=1;
#u(0):=0;
phi1[w]:=phi0;
N[bt]:=0.2;
mu1[w]:=mu(0);
k1[w]:=k(0);

               /                                     2\
eta -> mu1[bf] \1 + a[mu1] phi(eta) + b[mu1] phi(eta) /
              /                                   2\
eta -> k1[bf] \1 + a[k1] phi(eta) + b[k1] phi(eta) /
                              3880
                             998.2
                              773
                              4182
eta -> 2881.8 phi(eta) + 998.2
                    6                        6
       -1.1752324 10  phi(eta) + 4.1744724 10
eta -> ---------------------------------------
               2881.8 phi(eta) + 998.2        
              mu1[bf] (a[mu1] + 2 b[mu1] phi(eta))
                             39.11
                             533.9
                       0.000993000000000
                              7.47
                               0
                             0.597
                              0.5
                              0.04
                              0.2
           0.000993000000000 + 0.0388362300000 phi(0)

                                     2
              + 0.530162700000 phi(0)
                     0.597 + 4.45959 phi(0)
eq1:=subs(phi(0)=phi0,u(0)=0,eq1);
eq2:=subs(phi(0)=phi0,u(0)=0,eq2);
eq3:=subs(phi(0)=phi0,u(0)=0,eq3);
zet;

/  d   /  d         \\   (0.00339470952000)//                 
|----- |----- u(eta)|| +                    \0.000993000000000
\ deta \ deta       //                                        

                                                       2\   /
   + 0.0388362300000 phi(eta) + 0.530162700000 phi(eta) / + |
                                                            \

     1      /                                           /  d   
  ------- + |(0.0388362300000 + 1.06032540000 phi(eta)) |-----
  eta - 1   \                                           \ deta

          \\//                                            
  phi(eta)|| \0.000993000000000 + 0.0388362300000 phi(eta)
          //                                              

                            2\\ /  d         \
   + 0.530162700000 phi(eta) /| |----- u(eta)|
                              / \ deta       /
                                                  /          
                                                  |          
/  d   /  d         \\              1             |          
|----- |----- T(eta)|| + ------------------------ |0.7753836
\ deta \ deta       //   0.597 + 4.45959 phi(eta) \          

  /                                                              
  |              /             6                        6\       
  |2.66666666666 \-1.1752324 10  phi(eta) + 4.1744724 10 / u(eta)
  |--------------------------------------------------------------
  \                              p2                              

                   /  d           \
   + 5.75146288882 |----- phi(eta)|
                   \ deta         /

                                              /  d         \\\
     1.28968422855 (0.597 + 4.45959 phi(eta)) |----- T(eta)|||
                                              \ deta       /||
   - -------------------------------------------------------||
                             1 - eta                        //
    /  d           \                          /  d         \
    |----- phi(eta)| - 5.00000000000 phi(eta) |----- T(eta)|
    \ deta         /                          \ deta       /
                              0.5
res := dsolve({eq1=0,eq2=0,eq3=0,u(0)=0,D(u)(0)=a2,phi(0)=phi0,T(0)=0,D(T)(0)=1}, numeric,output=listprocedure,parameters=[a2,p2]):
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):
p:=proc(aa2,pp2) if not type([aa2,pp2],list(numeric)) then return 'procname(_passed)' end if:
res(parameters=[aa2,pp2]):
F0(1-zet);
end proc;
proc(aa2, pp2)  ...  end;


#for kkk from 0 to 2 do

papp2:=+4.364373285 : #+kkk*0.00001:

R0:=fsolve(p(AA,papp2*1e4)=0,AA=(0)..(1)):
print (R0);

R1:=evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet)):

R2:=papp2*1e4:
print (R1);
print (R2);
print (R1-R2);
print (seperate);
#end do:


                         0.190610828129
                   HFloat(43643.73285358669)
                          43643.73285
                  HFloat(3.586690581869334e-6)
                            seperate

ruu:=evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet)));;
phb:=evalf(2/(1-zet^2)*(Int((1-eta)*F1(eta),eta=0..1-zet)));
TTb:=evalf(2/(1-zet^2)*(Int((1-eta)*F2(eta),eta=0..1-zet)));
                   HFloat(0.0107135553557118)
                  HFloat(0.08288294476955627)
                   HFloat(0.1373974314161323)
with(plots):
res(parameters=[R0,R1]):
odeplot(res,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..zet);
odeplot(res,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..zet);




for giving the solution, I manually change the value of papp2 to obtaine R1-R2=0. I only ask you if it is possible, change this code in such a way that it automatically find the value of papp2. this code has no errors and converges.

thanks for your attention in advance.

Amir

Dear Carl

thanks for your reply

actually I can not understand your code appropraitely. (it is advanced)

But I Wrote a Code and found a solution to solve these equations.

But it doesnt find the _J_J_ automatically and I correct it manually.

the following code is my originla code

restart:
Digits := 12;
nta:=20;
gama1:=0:
phi0:=0.04;
                               12
                               20
                              0.04

rhocu:=2/(1-zet^2)*int((1-eta)*rho(eta)*c(eta)*u(eta),eta=0..1-zet);
#ff:=unapply((1-eta)*rho(eta)*c(eta)*u(eta),eta);
#rhocu:=2/(1-zet^2)*(1-zet)*( ff(0)/2 +sum(ff(0+ik*(1-zet)/nta),ik=1..nta-1))/nta;


eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])+((1/(eta-1)+1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
eq2:=diff(T(eta),eta,eta)+1/(k(eta)/k1[w])*(2/(1-zet^2)*rho(eta)*c(eta)*u(eta)/p2+( (a[k1]+2*b[k1]*phi(eta))/(1+a[k1]*phi1[w]+b[k1]*phi1[w]^2)*diff(phi(eta),eta)-k(eta)/k1[w]/(1-eta)*diff(T(eta),eta) ));
eq3:=diff(phi(eta),eta)-phi(eta)/(N[bt]*(1-gama1*T(eta))^2)*diff(T(eta),eta);
 2 (int((1 - eta) rho(eta) c(eta) u(eta), eta = 0 .. 1 - zet))
 -------------------------------------------------------------
                                  2                           
                           1 - zet                            
    /  d   /  d         \\   mu1[w]
    |----- |----- u(eta)|| + -------
    \ deta \ deta       //   mu(eta)

         /                 /  d           \\               
         |          mu_phi |----- phi(eta)||               
         |   1             \ deta         /| /  d         \
       + |------- + -----------------------| |----- u(eta)|
         \eta - 1           mu(eta)        / \ deta       /
                                /      /                        
                                |      |                        
/  d   /  d         \\     1    |      |2 rho(eta) c(eta) u(eta)
|----- |----- T(eta)|| + ------ |k1[w] |------------------------
\ deta \ deta       //   k(eta) |      |     /       2\         
                                \      \     \1 - zet / p2      

                                /  d           \
     (a[k1] + 2 b[k1] phi(eta)) |----- phi(eta)|
                                \ deta         /
   + -------------------------------------------
                                          2     
         1 + a[k1] phi1[w] + b[k1] phi1[w]      

            /  d         \\\
     k(eta) |----- T(eta)|||
            \ deta       /||
   - ---------------------||
        k1[w] (1 - eta)   ||
                          //
                                       /  d         \
                              phi(eta) |----- T(eta)|
           /  d           \            \ deta       /
           |----- phi(eta)| - -----------------------
           \ deta         /            N[bt]         
mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta);
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta);
rhop:=3880;
rhobf:=998.2;
cp:=773;
cbf:=4182;
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta);
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta);
mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta));

a[mu1]:=39.11;
b[mu1]:=533.9;
mu1[bf]:=9.93/10000;
a[k1]:=7.47;
b[k1]:=0;
k1[bf]:=0.597;
zet:=0.5;
#phi(0):=1;
#u(0):=0;
phi1[w]:=phi0;
N[bt]:=0.2;
mu1[w]:=mu(0);
k1[w]:=k(0);

               /                                     2\
eta -> mu1[bf] \1 + a[mu1] phi(eta) + b[mu1] phi(eta) /
              /                                   2\
eta -> k1[bf] \1 + a[k1] phi(eta) + b[k1] phi(eta) /
                              3880
                             998.2
                              773
                              4182
eta -> 2881.8 phi(eta) + 998.2
                    6                        6
       -1.1752324 10  phi(eta) + 4.1744724 10
eta -> ---------------------------------------
               2881.8 phi(eta) + 998.2        
              mu1[bf] (a[mu1] + 2 b[mu1] phi(eta))
                             39.11
                             533.9
                       0.000993000000000
                              7.47
                               0
                             0.597
                              0.5
                              0.04
                              0.2
           0.000993000000000 + 0.0388362300000 phi(0)

                                     2
              + 0.530162700000 phi(0)
                     0.597 + 4.45959 phi(0)
eq1:=subs(phi(0)=phi0,u(0)=0,eq1);
eq2:=subs(phi(0)=phi0,u(0)=0,eq2);
eq3:=subs(phi(0)=phi0,u(0)=0,eq3);
zet;

/  d   /  d         \\   (0.00339470952000)//                 
|----- |----- u(eta)|| +                    \0.000993000000000
\ deta \ deta       //                                        

                                                       2\   /
   + 0.0388362300000 phi(eta) + 0.530162700000 phi(eta) / + |
                                                            \

     1      /                                           /  d   
  ------- + |(0.0388362300000 + 1.06032540000 phi(eta)) |-----
  eta - 1   \                                           \ deta

          \\//                                            
  phi(eta)|| \0.000993000000000 + 0.0388362300000 phi(eta)
          //                                              

                            2\\ /  d         \
   + 0.530162700000 phi(eta) /| |----- u(eta)|
                              / \ deta       /
                                                  /          
                                                  |          
/  d   /  d         \\              1             |          
|----- |----- T(eta)|| + ------------------------ |0.7753836
\ deta \ deta       //   0.597 + 4.45959 phi(eta) \          

  /                                                              
  |              /             6                        6\       
  |2.66666666666 \-1.1752324 10  phi(eta) + 4.1744724 10 / u(eta)
  |--------------------------------------------------------------
  \                              p2                              

                   /  d           \
   + 5.75146288882 |----- phi(eta)|
                   \ deta         /

                                              /  d         \\\
     1.28968422855 (0.597 + 4.45959 phi(eta)) |----- T(eta)|||
                                              \ deta       /||
   - -------------------------------------------------------||
                             1 - eta                        //
    /  d           \                          /  d         \
    |----- phi(eta)| - 5.00000000000 phi(eta) |----- T(eta)|
    \ deta         /                          \ deta       /
                              0.5
res := dsolve({eq1=0,eq2=0,eq3=0,u(0)=0,D(u)(0)=a2,phi(0)=phi0,T(0)=0,D(T)(0)=1}, numeric,output=listprocedure,parameters=[a2,p2]):
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):
p:=proc(aa2,pp2) if not type([aa2,pp2],list(numeric)) then return 'procname(_passed)' end if:
res(parameters=[aa2,pp2]):
F0(1-zet);
end proc;
proc(aa2, pp2)  ...  end;


#for kkk from 0 to 2 do

papp2:=+4.364373285 : #+kkk*0.00001:

R0:=fsolve(p(AA,papp2*1e4)=0,AA=(0)..(1)):
print (R0);

R1:=evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet)):

R2:=papp2*1e4:
print (R1);
print (R2);
print (R1-R2);
print (seperate);
#end do:


                         0.190610828129
                   HFloat(43643.73285358669)
                          43643.73285
                  HFloat(3.586690581869334e-6)
                            seperate

ruu:=evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet)));;
phb:=evalf(2/(1-zet^2)*(Int((1-eta)*F1(eta),eta=0..1-zet)));
TTb:=evalf(2/(1-zet^2)*(Int((1-eta)*F2(eta),eta=0..1-zet)));
                   HFloat(0.0107135553557118)
                  HFloat(0.08288294476955627)
                   HFloat(0.1373974314161323)
with(plots):
res(parameters=[R0,R1]):
odeplot(res,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..zet);
odeplot(res,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..zet);




for giving the solution, I manually change the value of papp2 to obtaine R1-R2=0. I only ask you if it is possible, change this code in such a way that it automatically find the value of papp2. this code has no errors and converges.

thanks for your attention in advance.

Amir