In[]:=
clear[Re1,A1,Ha,λ,Pr,Rd,Ec,K1,Sc,S1,S2,β,pde1,pde2,pde3,pde4];
In[]:=
pde1[Re1_,A1_,Ha_,λ_]=D[u[t,x,y],t]+u[t,x,y]*(D[u[t,x,y],x])+v[t,x,y]*(D[u[t,x,y],y])D[D[u[t,x,y],y],y]+2*(A1^2/Re1)+((Ha+λ)*(A1-u[t,x,y]));​​pde2[λ_]=D[v[t,x,y],t]+u[t,x,y]*(D[v[t,x,y],x])+v[t,x,y]*(D[v[t,x,y],y])D[D[v[t,x,y],y],y]-λ*(v[t,x,y]);​​pde3[Re1_,Pr_,Rd_,Ec_,S1_]=D[theta[t,x,y],t]+u[t,x,y]*(D[theta[t,x,y],x])+v[t,x,y]*(D[theta[t,x,y],y])(1+(4/3)*Rd)/Pr*(D[D[theta[t,x,y],y],y])-(2/Re1)*(u[t,x,y])*(theta[t,x,y])-2*(S1/Re1)*(u[t,x,y])+Ec*((D[u[t,x,y],y])^2);​​pde4[Re1_,Sc_,K1_,S2_]=D[phi[t,x,y],t]+u[t,x,y]*(D[phi[t,x,y],x])+v[t,x,y]*(D[phi[t,x,y],y])(1/Sc)*(D[D[phi[t,x,y],y],y])-(2/Re1)*(u[t,x,y])*(phi[t,x,y])-(2/Re1)*S2*(u[t,x,y])-K1*(phi[t,x,y]);​​​​ics={u[0,x,y]0,v[0,x,y]0,theta[0,x,y]0,phi[0,x,y]0};​​With[{L=1},bcs[β_,A1_,S1_,S2_]={{u[t,0,y]0,v[t,0,y]0,theta[t,0,y]0,phi[t,0,y]0},{u[t,x,0]1,v[t,x,0]-β,theta[t,x,0]1-S1,phi[t,x,0]1-S2},u[t,x,L]A1,v[t,x,L]0,theta[t,x,L]0,phi[t,x,L]0}];
In[]:=
mol[n_Integer,o_:"Pseudospectral"]:={"MethodOfLines","SpatialDiscretization"{"TensorProductGrid","MaxPoints"n,"MinPoints"n,"DifferenceOrder"o}}​​mol[tf:False|True,sf_:Automatic]:={"MethodOfLines","DifferentiateBoundaryConditions"{tf,"ScaleFactor"sf}}
In[]:=
Clear@solfunc​​With[{pts=70,lb=1,L=1},solfunc[Re1_,Rd_,A1_:0.2,Ha_:1,tend_:2,λ_:2,Pr_:0.5,Ec_:1,K1_:0.2,Sc_:1,S1_:0.1,S2_:0.2,β_:0.3]:=NDSolveValue[{pde1[Re1,A1,Ha,λ],pde2[λ],pde3[Re1,Pr,Rd,Ec,S1],pde4[Re1,Sc,K1,S2],bcs[β,A1,S1,S2],ics},{u,v,theta,phi},{t,0,tend},{x,0,lb},{y,0,L},MethodUnion[mol[pts,4],mol[True,100]]]]​​
In[]:=
Table[sollst[Re1,Rd]=solfunc[Re1,Rd],{Re1,{0.2,0.4}},{Rd,{1,2,3}}];
NDSolveValue
:Warning: boundary and initial conditions are inconsistent.
NDSolveValue
:Warning: boundary and initial conditions are inconsistent.
NDSolveValue
:Warning: boundary and initial conditions are inconsistent.
General
:Further output of NDSolveValue::ibcinc will be suppressed during this calculation.
gr=Plot[Flatten@Table[-(x*(1+4/3*Rd)*(D[sollst[Re1,Rd][[3]][t,x,y],x]))-(x*(1+4/3*Rd)*(D[sollst[Re1,Rd][[3]][t,x,y],y]))/.y0,{Re1,{0.2,0.4}},{Rd,{1,2,3}}]//Evaluate,{t,0,1},PlotRangeAll,PlotStyle{Red,{Red,Thick,Dashed},Green,{Green,Thick,Dashed},Blue,{Blue,Thick,Dashed}},PlotLegendsPlaced[LineLegend[Automatic,{"Re1=0.2 ","Re1=0.25 "},LegendLayout"Column",LegendFunction(Framed[Style[#]]&),LegendMarkerSize{{30,1.3}}],{0.82,0.8}],FrameTrue,FrameStyleThick,FrameLabel->{Style["τ","Bold",14],Style["Nu","Bold",14]},LabelStyleBlack]
Out[]=
In[]:=
points=Cases[Normal@gr,Line[pts_,___]Flatten[pts,Depth[pts]-3],Infinity];
In[]:=
ifuncts=Interpolation[#,Method"Spline",InterpolationOrder2][x]&/@points;
In[]:=
data=Table[Prepend[ifuncts,x],{x,0,1,.01}];
In[]:=
TableForm[data,TableHeadings{None,Prepend[Array["f'",Length[points]],x]}];