## 195 Reputation

12 years, 133 days

Amir

## Convergence problem in my algorithm DSOL...

Hi,

I have the a code with some parameters including

Nr= 0, 50, 100

Ha=0, 5, 10

EPSILONE= 0, 0.5, 1

Phiavg= 0.02, 0.06, 0.1

0.1<NBT<10

I can give the solution for higher values of 5<NBT<10 and there is no problem. However, As I reduce the values of NBT, the convergence of the problem is hard. for some values of parameters I cannot find the solution. for example:

Nr=Ha=0

EPSILONE=1

Phiavg=0.06

NBT=0.3

I would be most grateful if you can tel me how change the algorithm to find the solution in the range of all parameters.

The code has been attached

code_7-8-2014_(1).mw

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?