> |
sys:=P*diff(f(eta), eta, eta, eta)/Q + lambda*theta(eta)*(S0 + S1*Qc*theta(eta)) + (n + 1)/2*(f(eta) + g(eta))*diff(f(eta), eta, eta) - n*diff(f(eta), eta)*(diff(f(eta), eta) + diff(g(eta), eta)) = 0, P*diff(g(eta), eta, eta, eta)/Q + (n + 1)/2*(f(eta) + g(eta))*diff(g(eta), eta, eta) - n*diff(g(eta), eta)*(diff(f(eta), eta) + diff(g(eta), eta)) = 0, (S2 - 2*R)*diff(theta(eta), eta, eta)/(S4*Pr) + (3*R)/(2*S4*Pr*(thetaw - 1))*diff((1 + (thetaw - 1)*theta(eta))^2, eta, eta) + (n + 1)/2*(f(eta) + g(eta))*diff(theta(eta), eta) = 0, f(0) = 0, g(0) = 0, theta(0) = 1, D(f)(0) = 1, D(g)(0) = c, D(f)(3) = 0, D(g)(3) = 0, theta(3) = 0;
|
> |
sys[3]; # Term 2 has a removable singularity at thetaw=1
|
> |
sys3:=map(normal,lhs(sys[3]))=0; # No singularity at thetaw=1
|
> |
isolate(sys3,diff(theta(eta),eta,eta));
|
When R > 0, theta(eta) > 0, and thetaw > 1 we see that the denominator is positive (in fact grater than R+S2).
> |
SYS:=sys[1..2],sys3,sys[4..-1];
|
> |
P := 0.7952865489: Q := 1.239340734: S4 := 1.007143462:S0 := 0.8330195568:S1 := 0.8330195568:
S2 := 1.087762421: Pr := 0.0157*2415/0.252: n:=3: lambda:=20: Qc:=2: c:=0.5:
|
> |
for s in SYS do s end do; # Viewing the system
|
Boundary value problems of this kind are notoriously difficult, so we begin with simply providing simple values for thetaw and R:
> |
res:=dsolve(eval({SYS},{thetaw=1.5,R=0.2}),numeric,initmesh=256,maxmesh=1024);
|
> |
dth0:=subs(res(0),diff(theta(eta),eta));
|
> |
plots:-odeplot(res,[[eta,f(eta)],[eta,g(eta)],[eta,theta(eta)]],0..3);
|
This procedure will be used for animation:
> |
resplot:=proc(thetaw0,R0) local res;
res:=dsolve(eval({SYS},{thetaw=thetaw0,R=R0}),numeric,_rest);
plots:-odeplot(res,[[eta,f(eta)],[eta,g(eta)],[eta,theta(eta)]],0..3,legend=[f,g,theta])
end proc;
|
> |
resplot(1.5,0.2,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]);
|
> |
plots:-animate(resplot,[1.5,R,initmesh=512,maxmesh=2048,
approxsoln=[f(eta)=0.4*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]],R=0.1..2,trace=24);
|
> |
plots:-animate(resplot,[thetaw,0.2,initmesh=512,maxmesh=2048,
approxsoln=[f(eta)=0.4*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta)]],thetaw=1..10,trace=24);
|
We notice that the general shape of the solutions doesn't change much with R and thetaw.
The following procedure takes values of thetaw and D(theta)(0) as input and outputs the corresponding value for R.
Notice the extra boundary condition D(theta)(0)=Dth0:
> |
Rproc:=proc(thetaw0,Dth0) local res;
if not [thetaw0,Dth0]::list(realcons) then return 'procname(_passed)' end if;
res:=dsolve(eval({SYS,D(theta)(0)=Dth0},thetaw=thetaw0),numeric,_rest);
subs(res(0),R)
end proc;
|
Here we should expect to see the result 0.2 for R even though the R value given in approxsoln is R = 1, and indeed we do:
> |
Rproc(1.5,dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]);
|
> |
Rproc(1.5,dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=2]);
|
> |
Rproc(1.5,-1,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]);
|
> |
plot(
Rproc(thetaw,-1,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]),
thetaw=1..2,labels=[thetaw,R]);
|
|