## 195 Reputation

9 years, 364 days

Amir

## Dslove numeric with know function...

Hi,

according to my previous question

http://www.mapleprimes.com/questions/201435-ODE-With-Constraint

I wrote the following code. at first, the code solve the equation for f and when it slves that, I want to solve Theta in such a way that use the values of f in previous calculation. I use the command 'known' but i couldnt find thesolution

I would be most grateful if you could help me in this problem

restart; # Notice that Restart (capital R) has no effect (to catch that use semicolon, not colon)
a:=0.13:
b:=0.41:
reynolds:=1.125*10^7;
Eq1:=diff(f(x),x\$3)+diff(f(x),x\$2)*f(x)+b^2*sqrt(2*reynolds)*diff(diff(f(x),x\$2)^2*x^2,x\$1);
Eq2:=diff(g(x),x\$3)+diff(g(x),x\$2)*g(x)+c*a^2*sqrt(2*reynolds)*diff(diff(g(x),x\$2)^2*x,x\$1);
eq1:=isolate(Eq1,diff(f(x),x,x,x));
eq2:=subs(g=f,isolate(Eq2,diff(g(x),x,x,x)));
EQ:=diff(f(x),x,x,x)=piecewise(x<c*0.1,rhs(eq1),rhs(eq2));

c:=75:
;
Q:=proc(pp2) local res,F0,F1,F2;
print(pp2);
if not type(pp2,numeric) then return 'procname(_passed)' end if:
res:=dsolve({EQ,f(0)=0,D(f)(0)=0,(D@@2)(f)(0)=pp2},numeric,output=listprocedure);
F0,F1,F2:=op(subs(subs(res),[f(x),diff(f(x),x),diff(f(x),x,x)])):
F1(c)-1;
end proc;

fsolve(Q(pp2)=0,pp2=(0..102));
se:=%;
res2:=dsolve({EQ,f(0)=0,D(f)(0)=0,(D@@2)(f)(0)=se},numeric,output=listprocedure);
G0,G1,G2:=op(subs(subs(res2),[f(x),diff(f(x),x),diff(f(x),x,x)])):

plots:-odeplot(res2,[seq([x,diff(f(x),[x\$i])],i=0..2)],0..2); #This plots from and past 0.1*c
pr:=1;
prt:=0.89;

Eq11:=diff(theta(x),x\$2)+pr*diff(theta(x),x\$1)*f(x)+pr/prt*b^2*sqrt(2*reynolds)*diff(diff(f(x),x\$2)*diff(theta(x),x\$1)*x^2,x\$1);
Eq22:=diff(g(x),x\$2)+pr*diff(g(x),x\$1)*f(x)+pr/prt*a^2*c*sqrt(2*reynolds)*diff(diff(f(x),x\$2)*diff(g(x),x\$1)*x^1,x\$1);
eq11:=isolate(Eq11,diff(f(x),x,x));
eq22:=subs(g=theta,isolate(Eq22,diff(g(x),x,x)));
EQT:=diff(theta(x),x,x)=piecewise(x<c*0.1,rhs(eq11),rhs(eq22));

QT:=proc(pp3) local res3,theta0,theta1;
print(pp3);
if not type(pp3,numeric) then return 'procname(_passed)' end if:
res3:=dsolve({EQT,theta(0)=1,D(theta)(0)=pp3,known=f},numeric,output=listprocedure);
theta0,theta1:=op(subs(subs(res),[theta(x),diff(theta(x),x)])):
theta0(c);
end proc;

fsolve(QT(pp3)=0,pp3=(0..200));
res3(0);

Amir

## ODE with constraint...

Hi

I have the following ODEs

Restart:
a:=0.13:
b:=0.41:
reynolds:=1.125*10^8:
Eq1:=diff(f(x),x\$3)+diff(f(x),x\$2)*f(x)+b^2*sqrt(2*reynolds)*diff(diff(f(x),x\$2)*f(x)*x^2,x\$1);
Eq2:=diff(g(x),x\$3)+diff(g(x),x\$2)*g(x)+c*a^2*sqrt(2*reynolds)*diff(diff(g(x),x\$2)*x,x\$1);
f(0)=0;
D(f(0)):=0;
# continuity condition
g(0.1*c)=f(0.1*c):
D(g(0.1*c))=D(f(0.1*c)):
(D@2)(g(0.1*c))=(D@2)(f(0.1*c)):

the value of c is unknown which must be obtained via D(g(c)):=1;

How can I solve it?

Amir

## widen the region of convergence ...

Dear collegues

I wrote the following code

restart:
Digits := 15;
a[k]:=0;
b[k]:=7.47;
a[mu]:=39.11;
b[mu]:=533.9;
mu[bf]:=9.93/10000;
k[bf]:=0.597;
ro[p]:=3880 ;
ro[bf]:= 998.2;
c[p]:= 773;
c[bf]:= 4182;
#mu[bf]:=1;
Gr[phi]:=0; Gr[T]:=0;
#dp:=0.1;
Ree:=1;
Pr:=1;
Nbt:=cc*NBTT+(1-cc^2)*6;

#######################
slip:=0.1;         ####
NBTT:=2;           ####
lambda:=0.1;       ####
phi_avg:=0.02;    ####
#######################

eq1:=diff( (1+a[mu]*phi(eta)+b[mu]*phi(eta)^2)*diff(u(eta),eta),eta)+dp/mu[bf]+Gr[T]*T(eta)-Gr[phi]*phi(eta);
eq2:=diff((1+a[k]*phi(eta)+b[k]*phi(eta)^2)*diff(T(eta),eta),eta)+lambda*T(eta)/k[bf];
eq3:=diff(phi(eta),eta)+1/Nbt*diff(T(eta),eta);
Q:=proc(pp2,fi0) local res,F0,F1,F2,a,INT0,INT10;
global Q1,Q2;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve({subs(dp=pp2,eq1)=0,eq2=0,eq3=0,u(0)=slip*D(u)(0),u(1)=-slip*D(u)(1),D(T)(0)=0,D(T)(1)=1,phi(0)=fi0}, numeric,output=listprocedure,continuation=cc);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):
INT0:=evalf(Int(F0(eta),eta=0..1));
INT10:=evalf(Int(F0(eta)*F1(eta),eta=0..1));
a[1]:=evalf(Int(F0(eta),eta=0..1))-Ree*Pr;;
a[2]:=INT10/INT0-phi_avg;
Q1(_passed):=a[1];
Q2(_passed):=a[2];
if type(procname,indexed) then a[op(procname)] else a[1],a[2] end if
end proc;
Q1:=proc(pp2,fi0) Q[1](_passed) end proc;
Q2:=proc(pp2,fi0) Q[2](_passed) end proc;
Optimization:-LSSolve([Q1,Q2],initialpoint=[0.3,0.0007]);

se:=%[2];
res2 := dsolve({subs(dp=se[1],eq1)=0,eq2=0,eq3=0,u(0)=slip*D(u)(0),u(1)=-slip*D(u)(1),D(T)(0)=0,D(T)(1)=1,phi(0)=se[2]}, numeric,output=listprocedure,continuation=cc);
G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):
TTb:=evalf(Int(G0(eta)*G2(eta)*(G1(eta)*ro[p]*c[p]+(1-G1(eta))*ro[bf]*c[bf] ),eta=0..1))/evalf(Int(G0(eta)*(G1(eta)*ro[p]*c[p]+(1-G1(eta))*ro[bf]*c[bf] ),eta=0..1));
with(plots):
odeplot(res2,[[eta,phi(eta)/phi_avg]],0..1);
odeplot(res2,[[eta,T(eta)/TTb]],0..1);
odeplot(res2,[[eta,u(eta)/(Ree*Pr)]],0..1);

res2(1);
Nuu:=(1/TTb);
1/((1+a[k]*G1(1)+b[k]*G1(1)^2)/(1+a[k]*phi_avg+b[k]*phi_avg^2));
(1/TTb)*(((1+a[k]*G1(1)+b[k]*G1(1)^2)/(1+a[k]*phi_avg+b[k]*phi_avg^2)));
>

I want to run the code for the value of NBTT in the range of 0.2 to 10. this code gave the results in the range of 4-10 easily. So, I used the continuation which improve the range of the results between 2-10. However, I coudnt gave the results when 0.2<NBTT<2. Would you please help me in this situation.

Also, It is to be said that the values of phi should be positive. in some ranges, I can see that phi(1) is negative. Can I place a condition in which the values phi restricted to be positive.

Amir

## Numerica Dsolve for Boundary values prob...

Maple 16

Hi,

I wrote the following code which is properly run

restart:

# parametrs

MUR:=(1-phi)^2.5:
RhoUR:=(1-phi+phi*rho[p]/rho[f]):
RhoCPR:=(1-phi+phi*rhocp[p]/rhocp[f]):
BetaUR:=(phi*rho[p]*beta[p]+(1-phi)*rho[f]*beta[f])/(RhoUR*rho[f])/beta[f]:

dqu3:=diff(h(x),x\$1)-RhoUR*BetaUR*T(x);
dqu2:=5*diff(T(x),x\$2)+k[f]/k[nf]*Pr*RhoCPR*f(x)*diff(T(x),x\$1);
dqu1:=5/(MUR)*diff(f(x),x\$3)
+ 2*(diff(h(x),x\$1)*x-h(x))
+RhoUR*(3*f(x)*diff(f(x),x\$2)-diff(f(x),x\$1)^2);
rho[f]:=998.2: cp[f]:=4182: k[f]:=0.597:   beta[f]:= 2.066/10000:
rho[p]:=3380: cp[p]:=773: k[p]:=36:   beta[p]:= 8.4/1000000:

k[nf]:=((k[p]+2*k[f])-2*phi*(k[f]-k[p]))/((k[p]+2*k[f])+phi*(k[f]-k[p])):
rhocp[nf]:=rho[p]*cp[p]*phi+rho[f]*cp[f]*(1-phi):
rhocp[p]:=rho[p]*cp[p]:
rhocp[f]:=rho[f]*cp[f]:

phi:=0.00:
binfinitive:=6: Pr:=7: lambda:=0:

with(plots):
pppe:=dsolve( {dqu1=0,dqu2=0,dqu3=0,T(0)=1,T(binfinitive)=0,f(0)=0,D(f)(0)=lambda,D(f)(binfinitive)=0,h(binfinitive)=0}, numeric );
-pppe(0);
print(odeplot(pppe,[x,diff(f(x),x)],0..binfinitive,color=black,numpoints=400));
print(odeplot(pppe,[[x,diff(f(x),x)]],0..binfinitive,color=black,numpoints=400));
print(odeplot(pppe,[[x,T(x)]],0..binfinitive,color=black,numpoints=400));

However, in some range of parameters, I must increase the value of binfinitive (for example binfinitive=50). however, my code is doesnt converge for higher values of 10 (at most). Can anyone change this algorithm in a way that it insensitive to the value of binfinitive?

Amir

## A Problem with LSSOLVE for ODE...

Maple 16

Following my previous question

http://www.mapleprimes.com/questions/200627-Lssolve-Midpoint

I wrote the following code

restart:
Phiavg:=0.06;
lambda:=0.05;
Ha:=0;
NBT:=0.5;
Nr:=500;
#N[bt]:=cc*NBT+(1-cc)*4; ## cc between 0 and 1
N[bt]:=cc*NBT+(1-cc^2)*0.75;

0.06
0.05
0
0.5
500
2
0.5 cc + 0.75 - 0.75 cc
eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])*(sigma-Nr*(phi(eta)-phi1[w])-(1-phi(eta))*T(eta)-Ha^2*u(eta))+((1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
eq2:=diff(T(eta),eta)-1/(k(eta)/k1[w]);
eq3:=diff(phi(eta),eta)-phi(eta)/(N[bt]*(1-gama1*T(eta))^2)*diff(T(eta),eta);
/  d   /  d         \\      1
|----- |----- u(eta)|| + ------- (mu1[w] (sigma - 500 phi(eta)
\ deta \ deta       //   mu(eta)

+ 500 phi1[w] - (1 - phi(eta)) T(eta)))

/  d           \ /  d         \
mu_phi |----- phi(eta)| |----- u(eta)|
\ deta         / \ deta       /
+ --------------------------------------
mu(eta)
/  d         \   k1[w]
|----- T(eta)| - ------
\ deta       /   k(eta)
/  d         \
phi(eta) |----- T(eta)|
/  d           \                       \ deta       /
|----- phi(eta)| - ----------------------------------------------
\ deta         /   /                       2\                   2
\0.5 cc + 0.75 - 0.75 cc / (1 - gama1 T(eta))
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)):
gama1:=0.00:
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:=1:
phi1[w]:=phi0:
mu1[w]:=mu(0):
k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,eq1);
eq2:=subs(phi(0)=phi0,eq2);
eq3:=subs(phi(0)=phi0,eq3);
/  d   /  d         \\   //
|----- |----- u(eta)|| + \\0.0009930000000 + 0.03883623000 phi0
\ deta \ deta       //

2\
+ 0.5301627000 phi0 / (sigma - 500 phi(eta) + 500 phi0

\//
- (1 - phi(eta)) T(eta))/ \0.0009930000000

2\
+ 0.03883623000 phi(eta) + 0.5301627000 phi(eta) / +

/                                       /  d           \ /  d
|(0.03883623000 + 1.060325400 phi(eta)) |----- phi(eta)| |-----
\                                       \ deta         / \ deta

\\//
u(eta)|| \0.0009930000000 + 0.03883623000 phi(eta)
//

2\
+ 0.5301627000 phi(eta) /
/  d         \     0.597 + 4.45959 phi0
|----- T(eta)| - ------------------------
\ deta       /   0.597 + 4.45959 phi(eta)
/  d         \
1. phi(eta) |----- T(eta)|
/  d           \               \ deta       /
|----- phi(eta)| - --------------------------
\ deta         /                           2
0.5 cc + 0.75 - 0.75 cc
Q:=proc(pp2,fi0) option remember; local res,F0,F1,F2,a,INT0,INT10,B;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if;
res := dsolve(subs(sigma=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}), numeric,output=listprocedure,initmesh=10, continuation=cc);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
INT0:=evalf(Int((abs(F0(eta)),eta=0..1)));
INT10:=evalf(Int(abs(F0(eta))*F1(eta),eta=0..1));
a[1]:=evalf(Int(F0(eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf),eta=0..1));
#a[1]:=evalf(Int((F0(eta),eta=0..1)));
a[2]:=(INT10/INT0-Phiavg)/Phiavg; #relative
[a[1],a[2]]
end proc:
Q1:=proc(pp2,fi0) Q(_passed)[1] end proc;
Q2:=proc(pp2,fi0) Q(_passed)[2] end proc;
proc(pp2, fi0)  ...  end;
proc(pp2, fi0)  ...  end;
#Q(116,0.0041);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[130,0.01]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[43.55,0.39]);
tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5.65,0.00036]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[12,0.003]); # khoob ba 1
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5,0.01]);
HFloat(5.65), HFloat(3.6e-4)
HFloat(5.650000070103341), HFloat(3.6e-4)
HFloat(5.65), HFloat(3.600105456508193e-4)
HFloat(29.63242379055208), HFloat(0.0205927592420527)
HFloat(12.803902258015825), HFloat(0.006395385884750864)
HFloat(12.803902403534572), HFloat(0.006395385884750864)
HFloat(12.803902258015825), HFloat(0.00639539649402585)
HFloat(12.804004931505949), HFloat(0.0063954867657199386)
HFloat(12.804107604996073), HFloat(0.006395587646689013)
HFloat(12.80400483062498), HFloat(0.006498160255844027)
HFloat(12.803902157134855), HFloat(0.006498059374874952)
HFloat(-1.0206939292143726), HFloat(-3.32764179807047e-4)
HFloat(-1.0206939079125088), HFloat(-3.32764179807047e-4)
HFloat(-1.0206939292143726), HFloat(-3.327536344433438e-4)
HFloat(18.749500943683863), HFloat(0.01993840615828979)
HFloat(3.9953780262640484), HFloat(0.00481041471606933)
HFloat(6.166152606930136), HFloat(0.00703619658484674)
HFloat(7.3193201827812295), HFloat(0.008218585352824569)
Error, (in Optimization:-LSSolve) complex value encountered
sigma:=tempe[2](1);
tempe[2](1)
phi0:=tempe[2](2);
tempe[2](2)
with(plots):

res2 := dsolve({eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}, numeric,output=listprocedure,continuation=cc);
Error, (in dsolve/numeric/process_input) boundary conditions specified at too many points: {0, 1, 2}, can only solve two-point boundary value problems
G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):
ruu:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet)))/(Phiavg*rhop+(1-Phiavg)*rhobf);
phb:=evalf((Int(abs(G0(eta))*G1(eta),eta=0..1))) / evalf((Int(abs(G0(eta)),eta=0..1))) ;
TTb:=evalf(Int(abs(G0(eta))*G2(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1))/evalf(Int(abs(G0(eta))*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1));
Error, invalid input: subs received res2, which is not valid for its 1st argument
/  /1.                                        \
| |                                           |
0.0008538922115 | |    |G0(eta)| (2881.8 G1(eta) + 998.2) deta|
| |                                           |
\/0.                                          /
/1.
|
|    |G0(eta)| G1(eta) deta
|
/0.
----------------------------
/1.
|
|
|    |G0(eta)| deta
/
0.
/Int(
1                               |
------------------------------------------------------------- |
/1.                                                         |
|                                                            \
|              /             6                       6\
|    |G0(eta)| \-1.1752324 10  G1(eta) + 4.1744724 10 / deta
/
0.

/             6                       6\ , eta = 0. .. 1.)
|G0(eta)| G2(eta) \-1.1752324 10  G1(eta) + 4.1744724 10 /

\
|
|
|
/
#rhouu:=evalf((Int((G1(eta)*rhop+(1-G1(eta))*rhobf)*G0(eta),eta=0..1)));

odeplot(res2,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..1);
#odeplot(res2,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..1);
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
rhou:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet))):

Nub:=(1/G2(1))*(((1+a[k1]*abs(G1(0))+b[k1]*abs(G1(0))^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2)));
0.6905123602 (1 + 7.47 |G1(0)|)
-------------------------------
G2(1)
(rhs(res2(0.0000000000001)[3])-rhs(res2(0)[3]))/0.0000000000001;
Error, invalid input: rhs received res2(0.1e-12)[3], which is not valid for its 1st argument, expr
sigma;
tempe[2](1)
NBT;
0.5
>

the above code has been worked for NBT=0.6 and higher, whereas as NBT decreases, the code doesnt converge easily.

How can I fix this problem?