## 195 Reputation

12 years, 134 days

Amir

Thanks Preben

Amir

## To Preben...

I get the correct results.

Thanks for your time and efforts for my code.

Amir

## a little change...

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

## response...

Wow!

It was very interesting. It is very quick.

Thanks a lot.

Amir

## response...

Wow!

It was very interesting. It is very quick.

Thanks a lot.

Amir

## BVP with constraining integrals...

I'm agree with you

## BVP with constraining integrals...

I'm agree with you

## New Question...

Hi

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?

## New Question...

Hi

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?

## @ carl...

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

## @ carl...

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

## the code...

Dear Carl

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.

Amir

## the code...

Dear Carl

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.