KIRAN SAJJAN

35 Reputation

4 Badges

0 years, 243 days

MaplePrimes Activity


These are replies submitted by KIRAN SAJJAN


 

restart:````

p1 := 0.1e-1; p2 := 0.2e-1; p3 := 0.1e-1; Px := p1+p2+p3

rf := 1050; kf := .52; cpf := 3617; sigmaf := .8

sigma1 := 25000; rs1 := 5200; ks1 := 6; cps1 := 670

sigma2 := 59.7*10^6; rs2 := 8933; ks2 := 400; cps2 := 385

sigma3 := 2380000; rs3 := 4250; ks3 := 8.9538; cps3 := 686.2

``

B1 := 1+2.5*Px+6.2*Px^2; B2 := 1+13.5*Px+904.4*Px^2; B3 := 1+37.1*Px+612.6*Px^2; B4 := (ks1+2*kf-2*Px*(kf-ks1))/(ks1+2*kf+Px*(kf-ks1)); B5 := (ks2+3.9*kf-3.9*Px*(kf-ks2))/(ks2+3.9*kf+Px*(kf-ks2)); B6 := (ks3+4.7*kf-4.7*Px*(kf-ks3))/(ks3+4.7*kf+Px*(kf-ks3))

a2 := B1*p1+B2*p2+B3*p3

a1 := 1-p1-p2-p3+p1*rs1/rf+p2*rs2/rf+p3*rs3/rf

a3 := 1-p1-p2-p3+p1*rs1*cps1/(rf*cpf)+p2*rs2*cps2/(rf*cpf)+p3*rs3*cps3/(rf*cpf)

a4 := B4*p1+B5*p2+B6*p3

``

a5 := 1+3*((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3)/(2+(p1*sigma1+p2*sigma2+p3*sigma3)/((p1+p2+p3)*sigmaf)-((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3))

NULL

NULL


ODEs:= [(a2+K)*(diff(U0(Y), Y$2))/a1-Ra*(diff(U0(Y), Y$1))+lambda0/a1-a5*M^2*U0(Y)/a1+K*(diff(N0(Y), Y$1))/a1+la*Ra*Theta0(Y)*(1+Qc*Theta0(Y)), (a2+K)*(diff(U1(Y), Y$2))/a1-H^2*l*U1(Y)-Ra*(diff(U1(Y), Y$1))+lambda1/a1-a5*M^2*U1(Y)/a1+K*(diff(N1(Y), Y$1))/a1+la*Ra*(Theta1(Y))(1+2*Qc*Theta0(Y)), diff(N0(Y), Y$2)-Ra*a1*Pj*(diff(N0(Y), Y$1))-2*n*N0(Y)-n*(diff(U0(Y), Y$1)), diff(N1(Y), Y$2)-Ra*a1*Pj*(diff(N1(Y), Y$1))-2*n*N1(Y)-n*(diff(U1(Y), Y$1))-H^2*a1*Pj*l*N1(Y), (a4/(a3*Pr)-delta*Ra^2/H^2+4*Rd*(1+(Tp-1)^3*Theta0(Y)^3+3*(Tp-1)^2*Theta0(Y)^2+(3*(Tp-1))*Theta0(Y))/(3*a3*Pr))*(diff(Theta0(Y), Y, Y))-Ra*(diff(Theta0(Y), Y))+a5*Ec*M^2*U0(Y)^2/a3+(a2+K)*Ec*(diff(U0(Y), Y))^2/a1+Q*Theta0(Y)/a3+4*(diff(Theta0(Y), Y))^2*Rd*(3*(Tp-1)+6*(Tp-1)^2*Theta0(Y)+3*(Tp-1)^3*Theta0(Y)^2)/(3*a3*Pr), (a4/(a3*Pr)-delta*Ra^2/H^2+4*Rd*(1+(Tp-1)^3*Theta0(Y)^3+3*(Tp-1)^2*Theta0(Y)^2+(3*(Tp-1))*Theta0(Y))/(3*a3*Pr))*(diff(Theta1(Y), Y, Y))-(H^2*l+2*Ra*l+Ra)*(diff(Theta1(Y), Y))+(Q/a3-delta*H^2*l^2)*Theta1(Y)+2*(a2+K)*Ec*(diff(U0(Y), Y))*(diff(U1(Y), Y))/a1+2*a5*Ec*M^2*U0(Y)*U1(Y)/a3+4*(diff(Theta0(Y), Y, Y))*Theta1(Y)*Rd*(3*(Tp-1)+6*(Tp-1)^2*Theta0(Y)+3*(Tp-1)^3*Theta0(Y)^2)/(3*a3*Pr)+4*Rd*(diff(Theta0(Y), Y))^2*(6*(Tp-1)^2*Theta1(Y)+6*(Tp-1)^3*Theta0(Y)*Theta1(Y))/(3*a3*Pr)+4*Rd*(diff(Theta1(Y), Y))*(diff(Theta0(Y), Y))*(6*(Tp-1)+6*(Tp-1)^3*Theta0(Y)^2+12*(Tp-1)^2*Theta0(Y))/(3*a3*Pr)
]:

``


(LB,UB):= (0,1):


BCs:= [
  
  U0(0) = 0, U1(0) = 0, N0(0) = 0, N1(0) = 0, Theta0(0) = 0, Theta1(0) = 0, U0(1) = 0, U1(1) = 0, N0(1) = 0, N1(1) = 0, Theta0(1) = 1, Theta1(1) = 0
]:

``


Params:= Record(
   
   M=  0.5, Rd=0.8,la=1.5,Tp=1.4,n=1.2,Q=0.4,Pj=0.001,Ra=0.4,Ec=0.5,    Pr= 21,   delta= 0.6,    t= (1/4)*Pi, lambda0=2,lambda1=3,   Qc= 0.1,    l= 1,K=0.4,H=3  ):
   

``


NBVs:= [   
   (a4+4*Rd*(1/3)+(4*Rd*(Tp-1)^3*(1/3))*(Theta0(0)+0.1e-2*exp(l*t)*Theta1(0)))*((D(Theta0))(0)+0.1e-2*exp(l*t)*(D(Theta1))(0)) = `Nu*`     # Nusselt numb
    
]:


Nu:= `Nu*`:


Solve:= module()
local
   nbvs_rhs:= rhs~(:-NBVs), #just the names
   Sol, #numeric dsolve BVP solution of any 'output' form
   ModuleApply:= subs(
      _Sys= {:-ODEs[], :-BCs[], :-NBVs[]},
      proc({
         M::realcons:=  Params:-M,
         Pr::realcons:= Params:-Pr,
         Rd::realcons:= Params:-Rd,
         la::realcons:= Params:-la,
         Tp::realcons:= Params:-Tp,
         n::realcons:= Params:-n,
         Q::realcons:= Params:-Q,
         Pj::realcons:= Params:-Pj,
         Ra::realcons:= Params:-Ra,
         Ec::realcons:= Params:-Ec,
         t::realcons:=  Params:-t,
         delta::realcons:= Params:-delta,
         lambda0::realcons:= Params:-lambda0,
         lambda1::realcons:= Params:-lambda1,
         Qc::realcons:= Params:-Qc,
         K::realcons:= Params:-K,
         l::realcons:= Params:-l,
         H::realcons:= Params:-H
         
      })
         Sol:= dsolve(_Sys, _rest, numeric);
         AccumData(Sol, {_options});
         Sol
      end proc
   ),
   AccumData:= proc(
      Sol::{Matrix, procedure, list({name, function}= procedure)},
      params::set(name= realcons)
   )
   local n, nbvs;
      if Sol::Matrix then
         nbvs:= seq(n = Sol[2,1][1,Pos(n)], n= nbvs_rhs)
      else
         nbvs:= (nbvs_rhs =~ eval(nbvs_rhs, Sol(:-LB)))[]
      fi;
      SavedData[params]:= Record[packed](params[], nbvs)
   end proc,
   ModuleLoad:= eval(Init);
export
   SavedData, #table of Records
   Pos, #Matrix column indices of nbvs
   Init:= proc()
      Pos:= proc(n::name) option remember; local p; member(n, Sol[1,1], 'p'); p end proc;
      SavedData:= table();
      return
   end proc ;
   ModuleLoad()
end module:

#procedure that generates 3-D plots (dropped-shadow contour + surface) of an expression


ParamPlot3d:= proc(
   Z::{procedure, `module`}, #procedure that extracts z-value from Solve's dsolve solution
   X::name= range(realcons), #x-axis-parameter range
   Y::name= range(realcons), #y-axis-parameter range
   FP::list(name= realcons), #fixed values of other parameters
   {
      #fraction of empty space above and below plot (larger "below"
      #value improves view of dropped-shadow contourplot):
      zmargin::[realcons,realcons]:= [.25,0.25],
      eta::realcons:= :-LB, #independent variable value
      dsolveopts::list({name, name= anything}):= [],
      contouropts::list({name, name= anything}):= [],
      surfaceopts::list({name, name= anything}):=[]    
   }
)
local
   LX:= lhs(X), RX:= rhs(X), LY:= lhs(Y), RY:= rhs(Y),
   Zremember:= proc(x,y)
   option remember; #Used because 'grid' should be the same for both plots.
      Z(
         Solve(
            LX= x, LY= y, FP[],
            #Default dsolve options can be changed by setting 'dsolveopts':
            'abserr'= 0.5e-20, 'interpolant'= false, 'output'= Array([eta]),  
            dsolveopts[]
         )
      )
   end proc,
   plotspec:= (Zremember, RX, RY),
   C:= plots:-contourplot(
      plotspec,
      #These default plot options can be changed by setting 'contouropts':
      'grid'= [25,25], 'contours'= 7, 'filled',
      'coloring'= ['yellow', 'orange'], 'color'= 'green',
      contouropts[]
   ),
   P:= plot3d(
      plotspec,
      #These default plot options can be changed by setting 'surfaceopts':
      'grid'= [25,25], 'style'= 'surfacecontour', 'contours'= 7,
      surfaceopts[]
   ),
   U, L #z-axis endpoints after margin adjustment
;
   #Stretch z-axis to include margins:
   (U,L):= ((Um,Lm,M,m)-> (M*(Lm-1)+m*Um, M*Lm+m*(Um-1)) /~ (Um+Lm-1))(
      zmargin[],
      (max,min)(op(3, indets(P, 'specfunc'('GRID'))[])) #actual z-axis range
   );
   plots:-display(
      [
         plots:-spacecurve(
            {
               [[lhs(RX),rhs(RY),U],[rhs(RX),rhs(RY),U],[rhs(RX),rhs(RY),L]], #yz backwall
               [[rhs(RX),rhs(RY),U],[rhs(RX),lhs(RY),U],[rhs(RX),lhs(RY),L]]  #xz backwall
            },
            'color'= 'grey', 'thickness'= 0
         ),
         plottools:-transform((x,y)-> [x,y,L])(C), #dropped-shadow contours
         P
      ],
      #These default plot options can be changed simply by putting the option in the
      #ParamPlot3d call:
      'view'= ['DEFAULT', 'DEFAULT', L..U], 'orientation'= [-135, 75], 'axes'= 'frame',
      'labels'= [lhs(X), lhs(Y), Z], 'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],
      'caption'= nprintf(cat("%a = %4.2f, "$nops(FP)-1, "%a = %4.2f"), (lhs,rhs)~(FP)[]),
      'captionfont'= ['TIMES', 14],
      'projection'= 2/3,   
      _rest
   )
end proc:

``

GetNu := proc (Sol::Matrix) options operator, arrow; Sol[2, 1][1, Solve:-Pos(:-Nu)] end proc

``

ParamPlot3d(
   GetNu,M = 0..1, K = 0..1, [      Pr= 21   ],
   labels= [M, K, Nu]);

Error, (in plot/iplot2d/levelcurve) could not evaluate expression

 

ParamPlot3d(GetNu, Ec = 0 .. 1, Q = 0 .. 2, [Pr = 21], labels = [Ec, Q, Nu])

Error, (in plot/iplot2d/levelcurve) could not evaluate expression

 

``

 

NULL


 

Download P6_3D_plots.mw

@tomleslie for the given equations implemented FDM method  i want to compare the graphs by dsolve and the FDM method for different values of Kp 

Actual plots are U versus Y  and T versus Y for different values of parameter kp =1,2,3

@tomleslie 

Mine actual equations I have given in above reply for that equations I want to solve by FDM and dsolve graphs for f versus  y and T versus y for the parameter value kp = 1,2,3

I know how draw by dsolve method but by FDM method I need to compare the plots

Please help me to solve 

In no post like this answeres posted 

@tomleslie for the given equations implemented FDM method  i want to compare the graphs by dsolve and the FDM method for different values of Kp 

Actual plots are U versus Y  and T versus Y for different values of parameter kp =1,2,3

Dear sir if i kept the boundary as 1 I m getting similar plots and table values for numerical and HPM but if i kept the boundary different like 3 or 4 plots are not converging By HPM method I have attached the code with boundary at 3 

Hpm_extension_paper.mw

@tomleslie 

I have changed the equations of theta_b and Q_b in the same loop but then also it is showing error for integration range Y=0..1

How to fix that error

Download error.mw

@tomleslie 

j is reffers to M values in the

Theta _b and Q_b

For same values I need to evaluate those Theta _b and Q_b values 

@mmcdara 

Thank you 

Plotting is correct sir but similar manner values of cf and nu for different values of parameters like M=1, cf=..... And Nu= .....

M=2, cf=.....and Nu=......

M=3,  cf=.....and Nu=......

 

@mmcdara  i have modified the terms 

in the plot i want different values of M=1, 2,3 With different colour and i want the table values for Cf and Nu for diffent values of M
 

``

restart; with(plots)

PDEtools[declare]((F, T, G, H)(Y), prime = Y)

` F`(Y)*`will now be displayed as`*F

 

` T`(Y)*`will now be displayed as`*T

 

` G`(Y)*`will now be displayed as`*G

 

` H`(Y)*`will now be displayed as`*H

 

`derivatives with respect to`*Y*`of functions of one variable will now be displayed with '`

(1)

p1 := 0.1e-1; p2 := 0.3e-1; p3 := 0.5e-1; p := p1+p2+p3

rf := 1050; kf := .52; cpf := 3617; sigmaf := .8

sigma1 := 25000; rs1 := 5200; ks1 := 6; cps1 := 670

sigma2 := 0.210e-5; rs2 := 5700; ks2 := 25; cps2 := 523

sigma3 := 6.30*10^7; rs3 := 10500; ks3 := 429; cps3 := 235

sigma4 := 10^(-10); rs4 := 3970; ks4 := 40; cps4 := 765

sigma5 := 1.69*10^7; rs5 := 7140; ks5 := 116; cps5 := 390

sigma6 := 4.10*10^7; rs6 := 19300; ks6 := 318; cps6 := 129

NULL

M := 1; S1 := .5; A := 1; delta := 0.1e-2; g := .1; Gr := .5; betu := .5; bett := .5; Pr := 21; Ec := .5; bet := 1; S2 := .5; Rd := 1; Q := .1; Ra := .5; S := .1

alp := 2

NULL

NULL

B1 := 1+2.5*p+6.2*p^2; B2 := 1+13.5*p+904.4*p^2; B3 := 1+37.1*p+612.6*p^2; B4 := (ks1+2*kf-2*p*(kf-ks1))/(ks1+2*kf+p*(kf-ks1)); B5 := (ks2+3.9*kf-3.9*p*(kf-ks2))/(ks2+3.9*kf+p*(kf-ks2)); B6 := (ks3+4.7*kf-4.7*p*(kf-ks3))/(ks3+4.7*kf+p*(kf-ks3)); B7 := (ks4+2*kf-2*p*(kf-ks4))/(ks4+2*kf+p*(kf-ks4)); B8 := (ks5+3.9*kf-3.9*p*(kf-ks5))/(ks5+3.9*kf+p*(kf-ks5)); B9 := (ks6+4.7*kf-4.7*p*(kf-ks6))/(ks6+4.7*kf+p*(kf-ks6))

a1 := B1*p1+B2*p2+B3*p3

a2 := 1-p1-p2-p3+p1*rs1/rf+p2*rs2/rf+p3*rs3/rf

a3 := 1-p1-p2-p3+p1*rs1*cps1/(rf*cpf)+p2*rs2*cps2/(rf*cpf)+p3*rs3*cps3/(rf*cpf)

a4 := B4*p1+B5*p2+B6*p3

NULL

a5 := 1+3*((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3)/(2+(p1*sigma1+p2*sigma2+p3*sigma3)/((p1+p2+p3)*sigmaf)-((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3))

NULL

NULL

a6 := B1*p1+B2*p2+B3*p3

a7 := 1-p1-p2-p3+p1*rs4/rf+p2*rs5/rf+p3*rs6/rf

a8 := 1-p1-p2-p3+p1*rs4*cps4/(rf*cpf)+p2*rs5*cps5/(rf*cpf)+p3*rs6*cps6/(rf*cpf)

a9 := B7*p1+B8*p2+B9*p3

NULL

a10 := 1+3*((p1*sigma4+p2*sigma5+p3*sigma6)/sigmaf-p1-p2-p3)/(2+(p1*sigma4+p2*sigma5+p3*sigma6)/((p1+p2+p3)*sigmaf)-((p1*sigma4+p2*sigma5+p3*sigma6)/sigmaf-p1-p2-p3))

W := sum(b[i]*Y^i, i = 0 .. 7); Theta := sum(c[i]*Y^i, i = 0 .. 7); U := sum(d[i]*Y^i, i = 0 .. 4); Phi := sum(h[i]*Y^i, i = 0 .. 4)

Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0]

 

Y^7*c[7]+Y^6*c[6]+Y^5*c[5]+Y^4*c[4]+Y^3*c[3]+Y^2*c[2]+Y*c[1]+c[0]

 

Y^4*d[4]+Y^3*d[3]+Y^2*d[2]+Y*d[1]+d[0]

 

Y^4*h[4]+Y^3*h[3]+Y^2*h[2]+Y*h[1]+h[0]

(2)

F := a1*(1+1/bet)*(diff(W, `$`(Y, 2)))+a2*Ra*(diff(W, Y))+A-a5*M*W-S2*W^2+a2*Gr*Theta-S*betu*(W-U) = 0

1-1.346703274*b[0]+.8111904760*b[1]+3.0560976*b[2]+.8111904760*c[1]*Y+.8111904760*c[2]*Y^2+.8111904760*c[3]*Y^3+.8111904760*c[4]*Y^4+.8111904760*c[5]*Y^5+.8111904760*c[6]*Y^6+.8111904760*c[7]*Y^7+64.1780496*Y^5*b[7]+45.8414640*Y^4*b[6]+30.5609760*Y^3*b[5]+18.3365856*Y^2*b[4]+9.1682928*Y*b[3]+5.678333332*Y^6*b[7]+4.867142856*Y^5*b[6]+4.055952380*Y^4*b[5]+3.244761904*Y^3*b[4]+2.433571428*Y^2*b[3]+1.622380952*Y*b[2]-1.346703274*b[1]*Y-1.346703274*b[2]*Y^2-1.346703274*b[3]*Y^3-1.346703274*b[4]*Y^4-1.346703274*b[5]*Y^5-1.346703274*b[6]*Y^6-1.346703274*b[7]*Y^7+0.5e-1*d[0]+0.5e-1*d[1]*Y+0.5e-1*d[2]*Y^2+0.5e-1*d[3]*Y^3+0.5e-1*d[4]*Y^4+.8111904760*c[0]-.5*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2 = 0

(3)

T := (a4+Rd)*(diff(Theta, `$`(Y, 2)))+a3*Pr*Ra*(diff(Theta, Y))+Q*Theta+Pr*alp*S*bett*(Theta-Phi)+Pr*Ec*((1+1/bet)*a1*(diff(W, Y))^2+a5*M*W^2+(1+1/bet)*a1*S1*W^2+S2*W^3+S*betu*(W-U)) = 0

.525*b[0]+2.20*c[1]*Y+2.20*c[2]*Y^2+2.20*c[3]*Y^3+2.20*c[4]*Y^4+2.20*c[5]*Y^5+2.20*c[6]*Y^6+2.20*c[7]*Y^7+.525*b[1]*Y+.525*b[2]*Y^2+.525*b[3]*Y^3+.525*b[4]*Y^4+.525*b[5]*Y^5+.525*b[6]*Y^6+.525*b[7]*Y^7-2.10*h[0]-.525*d[0]+47.59777865*Y^5*c[7]+33.99841332*Y^4*c[6]+22.66560888*Y^3*c[5]+13.59936533*Y^2*c[4]+6.799682664*Y*c[3]+71.67774537*Y^6*c[7]+61.43806746*Y^5*c[6]+51.19838955*Y^4*c[5]+40.95871164*Y^3*c[4]+30.71903373*Y^2*c[3]+20.47935582*Y*c[2]-2.10*Y^4*h[4]-2.10*Y^3*h[3]-2.10*Y^2*h[2]-2.10*Y*h[1]-.525*d[1]*Y-.525*d[2]*Y^2-.525*d[3]*Y^3-.525*d[4]*Y^4+2.266560888*c[2]+10.23967791*c[1]+2.20*c[0]+21.63764058*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+5.25*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^3+16.04451240*(7*Y^6*b[7]+6*Y^5*b[6]+5*Y^4*b[5]+4*Y^3*b[4]+3*Y^2*b[3]+2*Y*b[2]+b[1])^2 = 0

(4)

G := Ra*(diff(U, Y))+betu*(W-U) = 0

2.0*Y^3*d[4]+1.5*Y^2*d[3]+1.0*Y*d[2]+.5*d[1]+.5*b[7]*Y^7+.5*b[6]*Y^6+.5*b[5]*Y^5+.5*b[4]*Y^4-.5*d[4]*Y^4+.5*b[3]*Y^3-.5*d[3]*Y^3+.5*b[2]*Y^2-.5*d[2]*Y^2+.5*b[1]*Y-.5*d[1]*Y+.5*b[0]-.5*d[0] = 0

(5)

H := Ra*(diff(Phi, Y))+bett*(Theta-Phi) = 0

2.0*Y^3*h[4]+1.5*Y^2*h[3]+1.0*Y*h[2]+.5*h[1]+.5*c[7]*Y^7+.5*c[6]*Y^6+.5*c[5]*Y^5+.5*c[4]*Y^4-.5*Y^4*h[4]+.5*c[3]*Y^3-.5*Y^3*h[3]+.5*c[2]*Y^2-.5*Y^2*h[2]+.5*c[1]*Y-.5*Y*h[1]+.5*c[0]-.5*h[0] = 0

(6)

BCS := (D(W))(0) = 0, (D(Theta))(0) = 0, W(1) = -delta*(1+1/bet)*(D(W))(1), Theta(1) = 1+g*(D(Theta))(1), U(1) = -delta*(1+1/bet)*(D(W))(1), Phi(1) = 1+g*(D(Theta))(1)

NULL

W     := unapply(W, Y);Theta := unapply(Theta, Y);
U     := unapply(U, Y);Phi   := unapply(Phi, Y);
F     := unapply(F, Y);T     := unapply(T, Y);
G     := unapply(G, Y);H     := unapply(H, Y);



 

proc (Y) options operator, arrow; Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0] end proc

 

proc (Y) options operator, arrow; Y^7*c[7]+Y^6*c[6]+Y^5*c[5]+Y^4*c[4]+Y^3*c[3]+Y^2*c[2]+Y*c[1]+c[0] end proc

 

proc (Y) options operator, arrow; Y^4*d[4]+Y^3*d[3]+Y^2*d[2]+Y*d[1]+d[0] end proc

 

proc (Y) options operator, arrow; Y^4*h[4]+Y^3*h[3]+Y^2*h[2]+Y*h[1]+h[0] end proc

 

proc (Y) options operator, arrow; 1-1.346703274*b[0]+.8111904760*b[1]+3.0560976*b[2]+.8111904760*Y^7*c[7]+.8111904760*Y^6*c[6]+.8111904760*Y^5*c[5]+.8111904760*Y^4*c[4]+.8111904760*Y^3*c[3]+.8111904760*Y^2*c[2]+.8111904760*Y*c[1]+0.5e-1*Y^4*d[4]+0.5e-1*Y^3*d[3]+0.5e-1*Y^2*d[2]+0.5e-1*Y*d[1]-.5*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+64.1780496*Y^5*b[7]+45.8414640*Y^4*b[6]+30.5609760*Y^3*b[5]+18.3365856*Y^2*b[4]+9.1682928*Y*b[3]+5.678333332*Y^6*b[7]+4.867142856*Y^5*b[6]+4.055952380*Y^4*b[5]+3.244761904*Y^3*b[4]+2.433571428*Y^2*b[3]+1.622380952*Y*b[2]-1.346703274*Y^7*b[7]-1.346703274*Y^6*b[6]-1.346703274*Y^5*b[5]-1.346703274*Y^4*b[4]-1.346703274*Y^3*b[3]-1.346703274*Y^2*b[2]-1.346703274*Y*b[1]+0.5e-1*d[0]+.8111904760*c[0] = 0 end proc

 

proc (Y) options operator, arrow; .525*b[0]+2.20*Y^7*c[7]+2.20*Y^6*c[6]+2.20*Y^5*c[5]+2.20*Y^4*c[4]+2.20*Y^3*c[3]+2.20*Y^2*c[2]+2.20*Y*c[1]-.525*Y^4*d[4]-.525*Y^3*d[3]-.525*Y^2*d[2]-.525*Y*d[1]-2.10*Y^4*h[4]-2.10*Y^3*h[3]-2.10*Y^2*h[2]-2.10*Y*h[1]+21.63764058*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^2+.525*Y^7*b[7]+.525*Y^6*b[6]+.525*Y^5*b[5]+.525*Y^4*b[4]+.525*Y^3*b[3]+.525*Y^2*b[2]+.525*Y*b[1]+47.59777865*Y^5*c[7]+33.99841332*Y^4*c[6]+22.66560888*Y^3*c[5]+13.59936533*Y^2*c[4]+6.799682664*Y*c[3]+71.67774537*Y^6*c[7]+61.43806746*Y^5*c[6]+51.19838955*Y^4*c[5]+40.95871164*Y^3*c[4]+30.71903373*Y^2*c[3]+20.47935582*Y*c[2]+5.25*(Y^7*b[7]+Y^6*b[6]+Y^5*b[5]+Y^4*b[4]+Y^3*b[3]+Y^2*b[2]+Y*b[1]+b[0])^3+16.04451240*(7*Y^6*b[7]+6*Y^5*b[6]+5*Y^4*b[5]+4*Y^3*b[4]+3*Y^2*b[3]+2*Y*b[2]+b[1])^2-2.10*h[0]-.525*d[0]+2.266560888*c[2]+10.23967791*c[1]+2.20*c[0] = 0 end proc

 

proc (Y) options operator, arrow; 2.0*Y^3*d[4]+1.5*Y^2*d[3]+1.0*Y*d[2]+.5*d[1]+.5*Y^7*b[7]+.5*Y^6*b[6]+.5*Y^5*b[5]+.5*Y^4*b[4]-.5*Y^4*d[4]+.5*Y^3*b[3]-.5*Y^3*d[3]+.5*Y^2*b[2]-.5*Y^2*d[2]+.5*Y*b[1]-.5*Y*d[1]+.5*b[0]-.5*d[0] = 0 end proc

 

proc (Y) options operator, arrow; 2.0*Y^3*h[4]+1.5*Y^2*h[3]+1.0*Y*h[2]+.5*h[1]+.5*Y^7*c[7]+.5*Y^6*c[6]+.5*Y^5*c[5]+.5*Y^4*c[4]-.5*Y^4*h[4]+.5*Y^3*c[3]-.5*Y^3*h[3]+.5*Y^2*c[2]-.5*Y^2*h[2]+.5*Y*c[1]-.5*Y*h[1]+.5*c[0]-.5*h[0] = 0 end proc

(7)

 

z1 := (D(W))(0) = 0

z2 := W(1) = -delta*(1+1/bet)*(D(W))(1)

z3 := (D(Theta))(0) = 0

NULL

z4 := Theta(1) = 1+g*(D(Theta))(1)

 

z5 := U(1) = -delta*(1+1/bet)*(D(W))(1)

NULL

z6 := Phi(1) = 1+g*(D(Theta))(1)

z7 := F(0)

z8 := F(1)

NULL

z9 := (D(F))(0)

z10 := (D(F))(1)

z11 := (D(D(F)))(0)

z12 := (D(D(F)))(1)

NULL

z13 := T(0)

z14 := T(1)

z15 := (D(T))(0)

z16 := (D(T))(1)

z17 := (D(D(T)))(0)

z18 := (D(D(T)))(1)

z19 := G(0)

z20 := G(1)

z21 := H(0)

z22 := H(1)

z23 := (D(G))(0)

z24 := (D(G))(1)

z25 := (D(H))(0)

z26 := (D(H))(1)

NULL

``

``

NULL

Zs :=  [z1, z2, z3, z4, z5, z6, z7, z8, z9, z10, z11, z12, z13, z14,z15, z16,z17,z18,z19,z20,z21, z22, z23, z24,z25, z26]:
print~(Zs):

b[1] = 0

 

b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0] = -0.14e-1*b[7]-0.12e-1*b[6]-0.10e-1*b[5]-0.8e-2*b[4]-0.6e-2*b[3]-0.4e-2*b[2]-0.2e-2*b[1]

 

c[1] = 0

 

c[7]+c[6]+c[5]+c[4]+c[3]+c[2]+c[1]+c[0] = 1+.7*c[7]+.6*c[6]+.5*c[5]+.4*c[4]+.3*c[3]+.2*c[2]+.1*c[1]

 

d[4]+d[3]+d[2]+d[1]+d[0] = -0.14e-1*b[7]-0.12e-1*b[6]-0.10e-1*b[5]-0.8e-2*b[4]-0.6e-2*b[3]-0.4e-2*b[2]-0.2e-2*b[1]

 

h[4]+h[3]+h[2]+h[1]+h[0] = 1+.7*c[7]+.6*c[6]+.5*c[5]+.4*c[4]+.3*c[3]+.2*c[2]+.1*c[1]

 

1.-1.346703274*b[0]+.8111904760*b[1]+3.0560976*b[2]-.5*b[0]^2+0.5e-1*d[0]+.8111904760*c[0] = 0

 

1-1.346703274*b[0]-.5355127980*b[1]+3.331775278*b[2]+10.25516096*b[3]+20.23464423*b[4]+33.27022511*b[5]+49.36190359*b[6]+68.50967966*b[7]-.5*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2+0.5e-1*d[4]+0.5e-1*d[0]+0.5e-1*d[1]+0.5e-1*d[2]+0.5e-1*d[3]+.8111904760*c[7]+.8111904760*c[6]+.8111904760*c[5]+.8111904760*c[4]+.8111904760*c[3]+.8111904760*c[2]+.8111904760*c[1]+.8111904760*c[0] = 0

 

-1.346703274*b[1]+1.622380952*b[2]+9.1682928*b[3]+0.5e-1*d[1]+.8111904760*c[1]-1.0*b[0]*b[1] = 0

 

-1.346703274*b[1]-1.071025596*b[2]+9.995325838*b[3]+41.02064382*b[4]+101.1732212*b[5]+199.6213506*b[6]+345.5333251*b[7]+.20*d[4]+0.5e-1*d[1]+.10*d[2]+.15*d[3]+5.678333332*c[7]+4.867142856*c[6]+4.055952380*c[5]+3.244761904*c[4]+2.433571428*c[3]+1.622380952*c[2]+.8111904760*c[1]-1.0*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1]) = 0

 

-2.693406548*b[2]+4.867142856*b[3]+36.6731712*b[4]-1.0*b[1]^2-2.0*b[0]*b[2]+.10*d[2]+1.622380952*c[2] = 0

 

-2.693406548*b[2]-3.213076788*b[3]+39.98130333*b[4]+205.1032191*b[5]+607.0393269*b[6]+1397.349454*b[7]-1.0*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])+.60*d[4]+.10*d[2]+.30*d[3]+34.06999999*c[7]+24.33571428*c[6]+16.22380952*c[5]+9.734285712*c[4]+4.867142856*c[3]+1.622380952*c[2]-1.0*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2 = 0

 

.525*b[0]+21.63764058*b[0]^2+5.25*b[0]^3+16.04451240*b[1]^2-2.10*h[0]-.525*d[0]+2.266560888*c[2]+10.23967791*c[1]+2.20*c[0] = 0

 

.525*b[0]+.525*b[1]+.525*b[2]+.525*b[3]+.525*b[4]+.525*b[5]+.525*b[6]+.525*b[7]+5.25*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^3+21.63764058*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2-2.10*h[1]-2.10*h[2]-2.10*h[3]-2.10*h[4]-2.10*h[0]-.525*d[4]-.525*d[0]-.525*d[1]-.525*d[2]-.525*d[3]+121.4755240*c[7]+97.63648078*c[6]+76.06399843*c[5]+56.75807697*c[4]+39.71871639*c[3]+24.94591671*c[2]+12.43967791*c[1]+2.20*c[0]+16.04451240*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2 = 0

 

.525*b[1]+15.75*b[0]^2*b[1]+64.17804960*b[1]*b[2]-2.10*h[1]-.525*d[1]+6.799682664*c[3]+20.47935582*c[2]+2.20*c[1]+43.27528116*b[0]*b[1] = 0

 

.525*b[1]+1.050*b[2]+1.575*b[3]+2.100*b[4]+2.625*b[5]+3.150*b[6]+3.675*b[7]+15.75*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])+32.08902480*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])-2.10*h[1]-4.20*h[2]-6.30*h[3]-8.40*h[4]-2.100*d[4]-.525*d[1]-1.050*d[2]-1.575*d[3]+683.4553654*c[7]+456.3839906*c[6]+283.7903848*c[5]+158.8748656*c[4]+74.83775012*c[3]+24.87935582*c[2]+2.20*c[1]+43.27528116*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1]) = 0

 

1.050*b[2]+43.27528116*b[1]^2+86.55056232*b[0]*b[2]-4.20*h[2]-1.050*d[2]+27.19873066*c[4]+61.43806746*c[3]+4.40*c[2]+31.50*b[0]*b[1]^2+31.50*b[0]^2*b[2]+192.5341488*b[1]*b[3]+128.3560992*b[2]^2 = 0

 

1.050*b[2]+3.150*b[3]+6.300*b[4]+10.500*b[5]+15.750*b[6]+22.050*b[7]+43.27528116*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])-4.20*h[2]-12.60*h[3]-25.20*h[4]-6.300*d[4]-1.050*d[2]-3.150*d[3]+3194.687934*c[7]+1702.742309*c[6]+794.3743279*c[5]+299.3510005*c[4]+74.63806746*c[3]+4.40*c[2]+43.27528116*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2+31.50*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])^2+15.75*(b[7]+b[6]+b[5]+b[4]+b[3]+b[2]+b[1]+b[0])^2*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])+32.08902480*(7*b[7]+6*b[6]+5*b[5]+4*b[4]+3*b[3]+2*b[2]+b[1])*(210*b[7]+120*b[6]+60*b[5]+24*b[4]+6*b[3])+32.08902480*(42*b[7]+30*b[6]+20*b[5]+12*b[4]+6*b[3]+2*b[2])^2 = 0

 

.5*d[1]+.5*b[0]-.5*d[0] = 0

 

1.5*d[4]+1.0*d[3]+.5*d[2]+.5*b[7]+.5*b[6]+.5*b[5]+.5*b[4]+.5*b[3]+.5*b[2]+.5*b[1]+.5*b[0]-.5*d[0] = 0

 

.5*h[1]+.5*c[0]-.5*h[0] = 0

 

1.5*h[4]+1.0*h[3]+.5*h[2]+.5*c[7]+.5*c[6]+.5*c[5]+.5*c[4]+.5*c[3]+.5*c[2]+.5*c[1]+.5*c[0]-.5*h[0] = 0

 

1.0*d[2]+.5*b[1]-.5*d[1] = 0

 

4.0*d[4]+1.5*d[3]+3.5*b[7]+3.0*b[6]+2.5*b[5]+2.0*b[4]+1.5*b[3]+1.0*b[2]+.5*b[1]-.5*d[1] = 0

 

1.0*h[2]+.5*c[1]-.5*h[1] = 0

 

4.0*h[4]+1.5*h[3]+3.5*c[7]+3.0*c[6]+2.5*c[5]+2.0*c[4]+1.5*c[3]+1.0*c[2]+.5*c[1]-.5*h[1] = 0

(8)

IZs := indets(Zs);

numelems(Zs), numelems(IZs);

{b[0], b[1], b[2], b[3], b[4], b[5], b[6], b[7], c[0], c[1], c[2], c[3], c[4], c[5], c[6], c[7], d[0], d[1], d[2], d[3], d[4], h[0], h[1], h[2], h[3], h[4]}

 

26, 26

(9)

Zsol := fsolve(Zs)

{b[0] = .4528477729, b[1] = 0., b[2] = -.5456992939, b[3] = 0.9787892141e-1, b[4] = 0.5392589986e-1, b[5] = -.1004469472, b[6] = 0.5280749566e-1, b[7] = -0.9643937637e-2, c[0] = 1.688266265, c[1] = 0., c[2] = -2.715762562, c[3] = 8.023205978, c[4] = -18.19126346, c[5] = 22.70336194, c[6] = -13.86047071, c[7] = 3.251216630, d[0] = .2118383628, d[1] = -.2410094101, d[2] = -.1205047050, d[3] = .1233638337, d[4] = 0.2798182980e-1, h[0] = 1.242893264, h[1] = -.4453730008, h[2] = -.2226865004, h[3] = .4041352779, h[4] = -0.8041495805e-1}

(10)



F := unapply(eval(sum(b[i]*Y^i,i=0..7), Zsol), Y);
T := unapply(eval(sum(c[i]*Y^i,i=0..7), Zsol), Y);
G := unapply(eval(sum(d[i]*Y^i,i=0..4), Zsol), Y);
H := unapply(eval(sum(h[i]*Y^i,i=0..4), Zsol), Y);

proc (Y) options operator, arrow; -0.9643937637e-2*Y^7+0.5280749566e-1*Y^6-.1004469472*Y^5+0.5392589986e-1*Y^4+0.9787892141e-1*Y^3-.5456992939*Y^2+.4528477729 end proc

 

proc (Y) options operator, arrow; 3.251216630*Y^7-13.86047071*Y^6+22.70336194*Y^5-18.19126346*Y^4+8.023205978*Y^3-2.715762562*Y^2+1.688266265 end proc

 

proc (Y) options operator, arrow; 0.2798182980e-1*Y^4+.1233638337*Y^3-.1205047050*Y^2-.2410094101*Y+.2118383628 end proc

 

proc (Y) options operator, arrow; -0.8041495805e-1*Y^4+.4041352779*Y^3-.2226865004*Y^2-.4453730008*Y+1.242893264 end proc

(11)

``

plot(F(Y), Y = 0 .. 1, axes = boxed, color = green, thickness = 2, labels = [Y, W])

 




plot(T(Y), Y = 0 .. 1, axes = boxed, color = green, thickness = 2, labels = [Y, W])

 

plot(G(Y), Y = 0 .. 1, axes = boxed, color = green, thickness = 2, labels = [Y, W])

 

plot(H(Y), Y = 0 .. 1, axes = boxed, color = green, thickness = 2, labels = [Y, W])

 

Cf := a1*(1+1/bet)*(D(F))(1)

-1.275852826

(12)

Nu := (-a4-Rd)*(D(T))(1)

1.149666754

(13)

NULL

``


 

Download AGM.mw

@Carl Love

Dear sir i want 3d surface and contur plot 2d plots i m getting.

But according to your reference i kept the equation i m not able to get 3d plot for nusselt number 

@mmcdara 

without that can we plot 3d surface graph for Nusselt number Both nu and nu1 

BCs:= [(D(W))(0) = 0, (D(Theta))(0) = 0, W(1) = -delta*(1+1/bet)*(D(W))(1), Theta(1) = 1+g*(D(Theta))(1), V(1) = -delta*(1+1/bet)*(D(W))(1), Phi(1) = 1+g*(D(Theta))(1)
];
BCs1 := [(D(W1))(0) = 0, (D(Theta1))(0) = 0, W1(1) = -delta*(1+1/bet)*(D(W1))(1), Theta1(1) = 1+g*(D(Theta1))(1), V1(1) = -delta*(1+1/bet)*(D(W1))(1), Phi1(1) = 1+g*(D(Theta1))(1)]

BCs is for First Given Odes and BCs1 is for second system of odes

@mmcdara 

i have followed the below post according to that i have kept my equations

https://www.mapleprimes.com/posts/209715-Numerically-Solving-BVPs-That-Have-Many

@mmcdara  same problem i m using with different boundary conditios. it was exicuted in that i want 3D surface plot 

one error i got i am not able to rectify the error

@mmcdara  i want to plot 3d plots but getting some error

restart:

Digits:= trunc(evalhf(Digits)); #generally a very efficient setting

15

(1)

p1 := 0.1e-1; p2 := 0.3e-1; p3 := 0.5e-1; P := p1+p2+p3

rf := 1050; kf := .52; cpf := 3617; sigmaf := .8

sigma1 := 25000; rs1 := 5200; ks1 := 6; cps1 := 670

sigma2 := 0.210e-5; rs2 := 5700; ks2 := 25; cps2 := 523

sigma3 := 6.30*10^7; rs3 := 10500; ks3 := 429; cps3 := 235

sigma4 := 10^(-10); rs4 := 3970; ks4 := 40; cps4 := 765

sigma5 := 1.69*10^7; rs5 := 7140; ks5 := 116; cps5 := 390

sigma6 := 4.10*10^7; rs6 := 19300; ks6 := 318; cps6 := 129

NULL

B1 := 1+2.5*P+6.2*P^2; B2 := 1+13.5*P+904.4*P^2; B3 := 1+37.1*P+612.6*P^2; B4 := (ks1+2*kf-2*P*(kf-ks1))/(ks1+2*kf+P*(kf-ks1)); B5 := (ks2+3.9*kf-3.9*P*(kf-ks2))/(ks2+3.9*kf+P*(kf-ks2)); B6 := (ks3+4.7*kf-4.7*P*(kf-ks3))/(ks3+4.7*kf+P*(kf-ks3)); B7 := (ks4+2*kf-2*P*(kf-ks4))/(ks4+2*kf+P*(kf-ks4)); B8 := (ks5+3.9*kf-3.9*P*(kf-ks5))/(ks5+3.9*kf+P*(kf-ks5)); B9 := (ks6+4.7*kf-4.7*P*(kf-ks6))/(ks6+4.7*kf+P*(kf-ks6))

a1 := B1*p1+B2*p2+B3*p3

a2 := 1-p1-p2-p3+p1*rs1/rf+p2*rs2/rf+p3*rs3/rf

a3 := 1-p1-p2-p3+p1*rs1*cps1/(rf*cpf)+p2*rs2*cps2/(rf*cpf)+p3*rs3*cps3/(rf*cpf)

a4 := B4*p1+B5*p2+B6*p3

NULL

a5 := 1+3*((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3)/(2+(p1*sigma1+p2*sigma2+p3*sigma3)/((p1+p2+p3)*sigmaf)-((p1*sigma1+p2*sigma2+p3*sigma3)/sigmaf-p1-p2-p3))

NULL

NULL

a6 := B1*p1+B2*p2+B3*p3

a7 := 1-p1-p2-p3+p1*rs4/rf+p2*rs5/rf+p3*rs6/rf

a8 := 1-p1-p2-p3+p1*rs4*cps4/(rf*cpf)+p2*rs5*cps5/(rf*cpf)+p3*rs6*cps6/(rf*cpf)

a9 := B7*p1+B8*p2+B9*p3

NULL

a10 := 1+3*((p1*sigma4+p2*sigma5+p3*sigma6)/sigmaf-p1-p2-p3)/(2+(p1*sigma4+p2*sigma5+p3*sigma6)/((p1+p2+p3)*sigmaf)-((p1*sigma4+p2*sigma5+p3*sigma6)/sigmaf-p1-p2-p3))

NULL

 


ODEs:= [a1*(1+1/bet)*(diff(W(eta), eta, eta))+a2*Ra*(diff(W(eta), eta))+A-a5*M1*W(eta)-S2*W(eta)^2+a2*Gr*Theta(eta)-S*betu*(W(eta)-V(eta)) = 0, (a4+Rd)*(diff(Theta(eta), eta, eta))+a3*Pr*Ra*(diff(Theta(eta), eta))+Q*Theta(eta)+Pr*alp*S*bett*(Theta(eta)-Phi(eta))+Pr*Ec*((1+1/bet)*a1*(diff(W(eta), eta))^2+a5*M1*W(eta)^2+(1+1/bet)*a1*S1*W(eta)^2+S2*W(eta)^3+S*betu*(W(eta)-V(eta))) = 0, Ra*(diff(V(eta), eta))+betu*(W(eta)-V(eta)) = 0, Ra*(diff(Phi(eta), eta))+bett*(Theta(eta)-Phi(eta)) = 0
   
   
]:

<ODEs[]>:
   

ODEs1 := [a6*(1+1/bet)*(diff(W1(eta), eta, eta))+a7*Ra*(diff(W1(eta), eta))+A-a10*M1*W1(eta)-S2*W1(eta)^2+a7*Gr*Theta1(eta)-S*betu*(W1(eta)-V1(eta)) = 0, (a9+Rd)*(diff(Theta1(eta), eta, eta))+a8*Pr*Ra*(diff(Theta1(eta), eta))+Q*Theta1(eta)+Pr*alp*S*bett*(Theta1(eta)-Phi1(eta))+Pr*Ec*((1+1/bet)*a6*(diff(W1(eta), eta))^2+a10*M1*W1(eta)^2+(1+1/bet)*a6*S1*W1(eta)^2+S2*W1(eta)^3+S*betu*(W1(eta)-V1(eta))) = 0, Ra*(diff(V1(eta), eta))+betu*(W1(eta)-V1(eta)) = 0, Ra*(diff(Phi1(eta), eta))+bett*(Theta1(eta)-Phi1(eta)) = 0]; `<,>`(ODEs1[])

#lower and upper boundary values of the independent variable:
(LB,UB):= (0,1):

#boundary conditions:
BCs:= [(D(W))(0) = 0, (D(Theta))(0) = 0, W(1) = -delta*(1+1/bet)*(D(W))(1), Theta(1) = 1+g*(D(Theta))(1), V(1) = -delta*(1+1/bet)*(D(W))(1), Phi(1) = 1+g*(D(Theta))(1)
];

[(D(W))(0) = 0, (D(Theta))(0) = 0, W(1) = -delta*(1+1/bet)*(D(W))(1), Theta(1) = 1+g*(D(Theta))(1), V(1) = -delta*(1+1/bet)*(D(W))(1), Phi(1) = 1+g*(D(Theta))(1)]

(2)

BCs1 := [(D(W1))(0) = 0, (D(Theta1))(0) = 0, W1(1) = -delta*(1+1/bet)*(D(W1))(1), Theta1(1) = 1+g*(D(Theta1))(1), V1(1) = -delta*(1+1/bet)*(D(W1))(1), Phi1(1) = 1+g*(D(Theta1))(1)]

[(D(W1))(0) = 0, (D(Theta1))(0) = 0, W1(1) = -delta*(1+1/bet)*(D(W1))(1), Theta1(1) = 1+g*(D(Theta1))(1), V1(1) = -delta*(1+1/bet)*(D(W1))(1), Phi1(1) = 1+g*(D(Theta1))(1)]

(3)

 Params:= Record(
   
   S1=  0.5,   S2= 0.5, A=  1,
   Pr= 21, delta=0.01,g=0.1,
   Ec= 0.5, Ra=0.5,Q=1,Gr=0.5,Br=0.5,
   betu= 1.1,  bett= 1.1, S=  0.1,
   alp= 0.1,   
   bet= 2,
   M1= 1,
   Rd= 1   
):
   
   



NBVs:= [   
   a1*(1+1/bet)*D(W)(1) = `C*__f`,  # Skin friction coefficient
      -(a4+Rd)*D(Theta)(1) = `Nu*` ,     # Nusselt number
 
  
      a6*(1+1/bet)*D(W1)(1) = `C*__f1`,  # Skin friction coefficient
      -(a9+Rd)*D(Theta1)(1) = `Nu1*`      # Nusselt number
]:
Nu:= `Nu*`:
Cf:= `C*__f`:
Nu1:= `Nu1*`:
Cf1:= `C*__f1`:



Solve:= module()
local
   nbvs_rhs:= rhs~(:-NBVs), #just the names
   Sol, #numeric dsolve BVP solution of any 'output' form
   ModuleApply:= subs(
      _Sys= {:-ODEs[], :-BCs[], :-NBVs[]},
      proc({
         S::realcons:=  Params:-S,
         Pr::realcons:= Params:-Pr,
         Ec::realcons:= Params:-Ec,
         S1::realcons:= Params:-S1,
         M1::realcons:= Params:-M1,
         S2::realcons:= Params:-S2,
         A::realcons:= Params:-A,
         Rd::realcons:= Params:-Rd,
         g::realcons:=  Params:-g,
         delta::realcons:= Params:-delta,
         betu::realcons:= Params:-betu,
         Ra::realcons:= Params:-Ra,
         Q::realcons:= Params:-Q,
         Gr::realcons:= Params:-Gr,
         Br::realcons:= Params:-Br,
         alp::realcons:=  Params:-alp,
         bet::realcons:= Params:-bet,
         bett::realcons:= Params:-bett
      })
         Sol:= dsolve(_Sys, _rest, numeric);
         AccumData(Sol, {_options});
         Sol
      end proc
   ),
   AccumData:= proc(
      Sol::{Matrix, procedure, list({name, function}= procedure)},
      params::set(name= realcons)
   )
   local n, nbvs;
      if Sol::Matrix then
         nbvs:= seq(n = Sol[4,1][1,Pos(n)], n= nbvs_rhs)
      else
         nbvs:= (nbvs_rhs =~ eval(nbvs_rhs, Sol(:-LB)))[]
      fi;
      SavedData[params]:= Record[packed](params[], nbvs)
   end proc,
   ModuleLoad:= eval(Init);
export
   SavedData, #table of Records
   Pos, #Matrix column indices of nbvs
   Init:= proc()
      Pos:= proc(n::name) option remember; local p; member(n, Sol[1,1], 'p'); p end proc;
      SavedData:= table();
      return
   end proc ;
   ModuleLoad()
end module:
 



Solve:= module()
local
   nbvs_rhs:= rhs~(:-NBVs),
   Sol1,
   ModuleApply:= subs(
      _Sys1= {:-ODEs1[], :-BCs1[], :-NBVs[]},
      proc({
         S::realcons:=  Params:-S,
         Pr::realcons:= Params:-Pr,
         Ec::realcons:= Params:-Ec,
         S1::realcons:= Params:-S1,
         M1::realcons:= Params:-M1,
         S2::realcons:= Params:-S2,
         A::realcons:= Params:-A,
         Rd::realcons:= Params:-Rd,
         g::realcons:=  Params:-g,
         delta::realcons:= Params:-delta,
         betu::realcons:= Params:-betu,
         Ra::realcons:= Params:-Ra,
         Q::realcons:= Params:-Q,
         Gr::realcons:= Params:-Gr,
         Br::realcons:= Params:-Br,
         alp::realcons:=  Params:-alp,
         bet::realcons:= Params:-bet,
         bett::realcons:= Params:-bett
      })
         Sol1:= dsolve(_Sys1, _rest, numeric);
         AccumData(Sol1, {_options});
         Sol1
      end proc
   ),
   AccumData:= proc(
      Sol1::{Matrix, procedure, list({name, function}= procedure)},
      params::set(name= realcons)
   )
   local n, nbvs;
      if Sol1::Matrix then
         nbvs:= seq(n = Sol1[4,1][1,Pos(n)], n= nbvs_rhs)
      else
         nbvs:= (nbvs_rhs =~ eval(nbvs_rhs, Sol1(:-LB)))[]
      fi;
      SavedData[params]:= Record[packed](params[], nbvs)
   end proc,
   ModuleLoad:= eval(Init);
export
   SavedData,
   Pos1,
   Init:= proc()
      Pos1:= proc(n::name) option remember; local p; member(n, Sol1[1,1], 'p'); p end proc;
      SavedData:= table();
      return
   end proc ;
   ModuleLoad()
end module:

 

#procedure that generates 3-D plots (dropped-shadow contour + surface) of an expression


ParamPlot3d:= proc(
   Z::{procedure, `module`}, #procedure that extracts z-value from Solve's dsolve solution
   X::name= range(realcons), #x-axis-parameter range
   Y::name= range(realcons), #y-axis-parameter range
   FP::list(name= realcons), #fixed values of other parameters
   {
      #fraction of empty space above and below plot (larger "below"
      #value improves view of dropped-shadow contourplot):
      zmargin::[realcons,realcons]:= [.05,0.15],
      eta::realcons:= :-LB, #independent variable value
      dsolveopts::list({name, name= anything}):= [],
      contouropts::list({name, name= anything}):= [],
      surfaceopts::list({name, name= anything}):=[]    
   }
)
local
   LX:= lhs(X), RX:= rhs(X), LY:= lhs(Y), RY:= rhs(Y),
   Zremember:= proc(x,y)
   option remember; #Used because 'grid' should be the same for both plots.
      Z(
         Solve(
            LX= x, LY= y, FP[],
            #Default dsolve options can be changed by setting 'dsolveopts':
            'abserr'= 0.5e-7, 'interpolant'= false, 'output'= Array([eta]),  
            dsolveopts[]
         )
      )
   end proc,
   plotspec:= (Zremember, RX, RY),
   C:= plots:-contourplot(
      plotspec,
      #These default plot options can be changed by setting 'contouropts':
      'grid'= [25,25], 'contours'= 5, 'filled',
      'coloring'= ['yellow', 'orange'], 'color'= 'green',
      contouropts[]
   ),
   P:= plot3d(
      plotspec,
      #These default plot options can be changed by setting 'surfaceopts':
      'grid'= [25,25], 'style'= 'surfacecontour', 'contours'= 6,
      surfaceopts[]
   ),
   U, L #z-axis endpoints after margin adjustment
;
   #Stretch z-axis to include margins:
   (U,L):= ((Um,Lm,M,m)-> (M*(Lm-1)+m*Um, M*Lm+m*(Um-1)) /~ (Um+Lm-1))(
      zmargin[],
      (max,min)(op(3, indets(P, 'specfunc'('GRID'))[])) #actual z-axis range
   );
   plots:-display(  #final plot to display/return
      [
         plots:-spacecurve(
            {
               [[lhs(RX),rhs(RY),U],[rhs(RX),rhs(RY),U],[rhs(RX),rhs(RY),L]], #yz backwall
               [[rhs(RX),rhs(RY),U],[rhs(RX),lhs(RY),U],[rhs(RX),lhs(RY),L]]  #xz backwall
            },
            'color'= 'grey', 'thickness'= 0
         ),
         plottools:-transform((x,y)-> [x,y,L])(C), #dropped-shadow contours
         P
      ],
      #These default plot options can be changed simply by putting the option in the
      #ParamPlot3d call:
      'view'= ['DEFAULT', 'DEFAULT', L..U], 'orientation'= [-135, 75], 'axes'= 'frame',
      'labels'= [lhs(X), lhs(Y), Z], 'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],
      'caption'= nprintf(cat("%a = %4.2f, "$nops(FP)-1, "%a = %4.2f"), (lhs,rhs)~(FP)[]),
      'captionfont'= ['TIMES', 14],
      'projection'= 2/3,   
      _rest
   )
end proc:

GetNu := proc (Sol::Matrix) options operator, arrow; Sol[4, 1][1, Solve:-Pos(:-Nu)] end proc; GetNu1 := proc (Sol1::Matrix) options operator, arrow; Sol1[4, 1][1, Solve1:-Pos(:-Nu1)] end proc

 

ParamPlot3d(
   GetNu,S= 0..5, M1= 0..5, [S1=  0.5,   S2= 0.5, A=  1,
   Pr= 21, delta=0.01,g=0.1,
   Ec= 0.5, Ra=0.5,Q=1,Gr=0.5,Br=0.5,
   betu= 1.1,  bett= 1.1,
   alp= 0.1,   
   bet= 2,
   
   Rd= 1   ],
   labels= [S, M1, Nu]
);

Error, (in plot/iplot2d/levelcurve) could not evaluate expression

 

ParamPlot3d(GetNu1, S = 0 .. 5, M1 = 0 .. 5, [S1 = .5, S2 = .5, A = 1, Pr = 21, delta = 0.1e-1, g = .1, Ec = .5, Ra = .5, Q = 1, Gr = .5, Br = .5, betu = 1.1, bett = 1.1, alp = .1, bet = 2, Rd = 1], labels = [S, M1, Nu])

Error, (in plot/iplot2d/levelcurve) could not evaluate expression

 

NULL

Download 3d_plots.mw

restart; with(plots)

_local(gamma); _local(Re)

p1 := 0.1e-1; p2 := 0.3e-1; p3 := 0.5e-1; p := p1+p2+p3

rf := 1050; kf := .52; cpf := 3617; `&sigma;f` := .8

sigma1 := 25000; rs1 := 5200; ks1 := 6; cps1 := 670

sigma2 := 2*10^(-6); rs2 := 5700; ks2 := 25; cps2 := 523

sigma3 := 6.30*10^7; rs3 := 10500; ks3 := 429; cps3 := 235

``

S1 := 1; S2 := 1; A := 1; Pr := 21; delta := 1; Rd := 1; gamma := 1; Ec := 1; Re := 1; Q := 1; Gr := 1; Br := 1; `&beta;u` := 1; `&beta;t` := 1; S := 1; alpha := 1; beta := 1

NULL

``

``

B1 := 1+2.5*p+6.2*p^2; B2 := 1+13.5*p+904.4*p^2; B3 := 1+37.1*p+612.6*p^2; B4 := (ks1+2*kf-2*p*(kf-ks1))/(ks1+2*kf+p*(kf-ks1)); B5 := (ks2+3.9*kf-3.9*p*(kf-ks2))/(ks2+3.9*kf+p*(kf-ks2)); B6 := (ks3+4.7*kf-4.7*p*(kf-ks3))/(ks3+4.7*kf+p*(kf-ks3))

a1 := B1*p1+B2*p2+B3*p3

a2 := 1-p1-p2-p3+p1*rs1/rf+p2*rs2/rf+p3*rs3/rf

a3 := 1-p1-p2-p3+p1*rs1*cps1/(rf*cpf)+p2*rs2*cps2/(rf*cpf)+p3*rs3*cps3/(rf*cpf)

a4 := B4*p1+B5*p2+B6*p3

NULL

a5 := 1+3*((p1*sigma1+p2*sigma2+p3*sigma3)/`&sigma;f`-p1-p2-p3)/(2+(p1*sigma1+p2*sigma2+p3*sigma3)/((p1+p2+p3)*`&sigma;f`)-((p1*sigma1+p2*sigma2+p3*sigma3)/`&sigma;f`-p1-p2-p3))

``

NULL

OdeSys := a1*(1+1/beta)*(diff(W(eta), eta, eta))+a2*Re*(diff(W(eta), eta))+A-a5*M*W(eta)-S2*W(eta)^2+a2*Gr*Theta(eta)-S*`&beta;u`*(W(eta)-Wd(eta)), (a4+Rd)*(diff(Theta(eta), eta, eta))+a3*Pr*Re*(diff(Theta(eta), eta))+Q*Theta(eta)+Pr*alpha*S*`&beta;t`*(Theta(eta)-Theta_d(eta))+Pr*Ec*{(1+1/beta)*a1*(diff(W(eta), eta))^2+a5*M*W(eta)^2+(1+1/beta)*a1*S1*W(eta)^2+S2*W(eta)^3+S*`&beta;u`*(W(eta)-Wd(eta))}, Re*(diff(Wd(eta), eta))+`&beta;u`*(W(eta)-Wd(eta)), Re*(diff(Theta_d(eta), eta))+`&beta;t`*(Theta(eta)-Theta_d(eta)); Cond := (D(W))(0) = 0, (D(Theta))(0) = 0, (D(Wd))(0) = 0, (D(Theta_d))(0) = 0, W(1) = -delta*(1+1/beta)*(D(W))(1), Theta(1) = 1+gamma*(D(Theta))(1), Wd(1) = -delta*(1+1/beta)*(D(Wd))(1), Theta_d(1) = 1+gamma*(D(Theta_d))(1)

1.5280488*(diff(diff(W(eta), eta), eta))+1.622380952*(diff(W(eta), eta))+1-1.296703274*M*W(eta)-W(eta)^2+1.622380952*Theta(eta)-W(eta)+Wd(eta), 1.133280444*(diff(diff(Theta(eta), eta), eta))+20.47935582*(diff(Theta(eta), eta))+22*Theta(eta)-21*Theta_d(eta)+21*{1.5280488*(diff(W(eta), eta))^2+1.296703274*M*W(eta)^2+1.5280488*W(eta)^2+W(eta)^3+W(eta)-Wd(eta)}, diff(Wd(eta), eta)+W(eta)-Wd(eta), diff(Theta_d(eta), eta)+Theta(eta)-Theta_d(eta)

 

(D(W))(0) = 0, (D(Theta))(0) = 0, (D(Wd))(0) = 0, (D(Theta_d))(0) = 0, W(1) = -2*(D(W))(1), Theta(1) = 1+(D(Theta))(1), Wd(1) = -2*(D(Wd))(1), Theta_d(1) = 1+(D(Theta_d))(1)

(1)

``

Error, (in dsolve/numeric/bvp/convertsys) too many boundary conditions: expected 6, got 8

 

MVals := [1, 3, 5]

for j to numelems(MVals) do Ans[j] := dsolve(eval([OdeSys, Cond], M = MVals[j]), numeric, output = listprocedure) end do

Error, (in dsolve/numeric/bvp/convertsys) too many boundary conditions: expected 6, got 8

 

 

``

``

interface(rtablesize = 100); interface(displayprecision = 5); Matrix([[0, C1, C1, C1, C2, C2, C2, C3, C3, C3], [eta, Cf, Nu, Ng, Cf, Nu, Ng, Cf, Nu, Ng], seq([j, seq([a1*(1+1/beta)*(eval(diff(W(eta), eta), Ans[k]))(j), -(a4+Rd)*(eval(diff(Theta(eta), eta), Ans[k]))(j), (a4+Rd)*(eval((diff(Theta(eta), eta))^2, Ans[k]))(j)+a5*M*Br*(eval(W(eta)^2, Ans[k]))(j)+Br*(1+1/beta)*a1*(eval((diff(W(eta), eta))^2, Ans[k]))(j)+a1*S1*Br*(1+1/beta)^(eval(W(eta)^2, Ans[k]))(j)][], k = 1 .. numelems(MVals))], j = 1)]); interface(rtablesize = 10); interface(displayprecision = -1)

Error, invalid input: numelems expects its 1st argument, t, to be of type indexable, but received RdVals

 

NULL

with(plots):
  cols := [red, blue, black]:

 plotA:= display
  ( [ seq
      ( odeplot
        ( Ans[k],[eta,D(W(eta))],
          eta=0..1,
          color=cols[k]
        ),
        k=1..numelems(MVals)
      )
    ],
    'axes'= 'boxed',labels=[eta,'W(eta)'],title= Typesetting:-mtext(" Velocity distribution" ,color= Blue)
  );
 

with(plots):
  cols := [red, blue, black]:

plotB:= display( [ seq( odeplot
        ( Ans[k],[eta,Theta(eta)],
          eta=0..1,
          color=cols[k]
        ),
        k=1..numelems(MVals)
      )
    ],
    'axes'= 'boxed',labels=[eta,'Theta(eta)'],title= Typesetting:-mtext("Temperature distribution" ,color= Blue)
  );

Error, invalid input: numelems expects its 1st argument, t, to be of type indexable, but received RdVals

 

Error, invalid input: numelems expects its 1st argument, t, to be of type indexable, but received RdVals

 

 with(plots):
  cols := [red, blue, black]:

 plotC:= display
  ( [ seq
      ( odeplot
        ( Ans[k],[eta,D(Wd(eta))],
          eta=0..1,
          color=cols[k]
        ),
        k=1..numelems(MVals)
      )
    ],
   'axes'= 'boxed','linestyle' = 'dashdot',labels=[eta,'Wd(eta)'],title= Typesetting:-mtext("Velocity distribution" ,color= Red)
  );

Error, invalid input: numelems expects its 1st argument, t, to be of type indexable, but received RdVals

 

with(plots); cols := [red, blue, black]; plotD := display*([seq*(odeplot*(Ans[k], [eta, Theta_d(eta)], eta = 0 .. 1, color = cols[k]), k = 1 .. numelems(MVals))], 'axes' = 'boxed', 'linestyle' = 'dashdot', labels = [eta, 'Theta_d(eta)'], title = Typesetting:-mtext("Velocity distribution", color = Red))

Error, invalid input: numelems expects its 1st argument, t, to be of type indexable, but received RdVals

 

plots:-display([plotA, plotC])

Error, (in plots:-display) expecting plot structures but received: [plotA, plotC]

 

plots:-display([plotB, plotD])

Error, (in plots:-display) expecting plot structures but received: [plotB, plotD]

 

````

Download dust_particle_shape_hybrid_NF.mw

1 2 3 Page 1 of 3