KIRAN SAJJAN

followed by below given post...

Answered: KIRAN SAJJAN 40

May 15 2023

0 0

@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

yes...

Answered: KIRAN SAJJAN 40

May 15 2023

0 0

@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

I have changed the boundary...

Answered: KIRAN SAJJAN 40

May 15 2023

0 0

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

>

restart:

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

(1)

>

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))

>

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))

>

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[]>:

>

>	#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

Download 3d_plots.mw

How to fix this error...

Answered: KIRAN SAJJAN 40

May 11 2023

0 0

>

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))

>

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

>

$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*`β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*`β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*`βu`*(W(eta)-Wd(eta))}, Re*(diff(Wd(eta), eta))+`βu`*(W(eta)-Wd(eta)), Re*(diff(Theta_d(eta), eta))+`β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)$

(1)

>

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

>

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

>

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

>

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

>

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

>

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

>

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

Download dust_particle_shape_hybrid_NF.mw

now executing sir...

Answered: KIRAN SAJJAN 40

May 05 2023

0 0

@tomleslie now i want similar graph from desolve and table values

i dont know the command to plot the graph i want to compair the results Download FDM_practice.mw

>

restart; with(plots)

>

a1 := 1/(1-p1-p2)^2.5

>

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

>

a4 := ((ks1*p1+ks2*p2)/(p1+p2)+2*kf+2*(ks1*p1+ks2*p2)-2*(p1+p2)*kf)/((ks1*p1+ks2*p2)/(p1+p2)+2*kf-ks1*p1-ks2*p2+(p1+p2)*kf)

>

a5 := 1+3*((p1*sigma1+p2*sigma2)/`σf`-p1-p2)/(2+(p1*sigma1+p2*sigma2)/((p1+p2)*`σf`)-((p1*sigma1+p2*sigma2)/`σf`-p1-p2))

>

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

>

interface(rtablesize = 100); interface(displayprecision = 5); Matrix([[Y, Nu1a, `Ν1b`, Nu1c], seq([j, seq([(a4+4*Rd*(1/3))*((eval(diff(Theta0(Y), Y), Ans[k]))(j)+0.1e-2*evalf(exp(l*t))*(eval(diff(Theta1(Y), Y), Ans[k]))(j))][], k = 1 .. numelems(RVals))], j = 0)]); interface(rtablesize = 10); interface(displayprecision = -1)

Matrix(%id = 18446746568164137670)

(1)

>

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

plotA:= display
  ( [ seq
      ( odeplot
        ( Ans[k],[Y,U0(Y)],
          Y=0..1,
          color=cols[k]
        ),
        k=1..numelems(RVals)
      )
    ],
    'axes'= 'boxed',labels=[Y,'U[0]'],title= Typesetting:-mtext("Steady Velocity distribution" ,color= Red)
  );
plotB:= display
  ( [ seq
      ( odeplot
        ( Ans[k],[Y,U1(Y)],
          Y=0..1,
          color=cols[k]
        ),
        k=1..numelems(RVals)
      )
    ],
    'axes'= 'boxed',labels=[Y,'U[1]'],title= Typesetting:-mtext("Unsteady Velocity distribution" ,color= Blue)

  );

>

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

plotE:= display( [ seq( odeplot
        ( Ans[k],[Y,Theta0(Y)],
          Y=0..1,
          color=cols[k]
        ),
        k=1..numelems(RVals)
      )
    ],
    'axes'= 'boxed',labels=[Y,'Theta[0]'],title= Typesetting:-mtext("Steady Temperature distribution" ,color= Red)
  );
plotF:= display
  ( [ seq
      ( odeplot
        ( Ans[k],[Y,Theta1(Y)],
          Y=0..1,
          color=cols[k]
        ),
        k=1..numelems(RVals)
      )
    ],
    'axes'= 'boxed',labels=[Y,'Theta[1]'],title= Typesetting:-mtext("Unsteady Temperature distribution" ,color= Blue)

  );

>

(2)

>

restart;

>

a1 := 1/(1-p1-p2)^2.5

>

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

>

a4 := ((ks1*p1+ks2*p2)/(p1+p2)+2*kf+2*(ks1*p1+ks2*p2)-2*(p1+p2)*kf)/((ks1*p1+ks2*p2)/(p1+p2)+2*kf-ks1*p1-ks2*p2+(p1+p2)*kf)

>

a5 := 1+3*((p1*sigma1+p2*sigma2)/`σf`-p1-p2)/(2+(p1*sigma1+p2*sigma2)/((p1+p2)*`σf`)-((p1*sigma1+p2*sigma2)/`σf`-p1-p2))

>

a := 0; b := 1; N := 9; h := (b-a)/(N+1)

>

# Initial conditions
for i from 0 to N do
  U[i] := h*i+1:
end do:

for i from 0 to N do
  T[i] := h*i+1:
end do:
for i from 0 to N do
  C[i] := h*i+1:
end do:

# Boundary conditions

  U0[0] := 0:
  U0[N+1] := 0:
  U1[0] := 0:
  U1[N+1] := 0:
  Theta0[0] := 0:
  Theta0[N+1] := 1:
  Theta1[0] := 0:
  Theta1[N+1] := 0:

d1:= (U0[i+1]-U0[i-1])/(2*h):
d2:=(U0[i+1]-2*U0[i]+U0[i-1])/h^2:
d3:= (U1[i+1]-U1[i-1])/(2*h):d4:=(U1[i+1]-2*U1[i]+U1[i-1])/h^2:
d5:= (Theta0[i+1]-Theta0[i-1])/(2*h):d6:=(Theta0[i+1]-2*Theta0[i]+Theta0[i-1])/h^2:
d7:= (Theta1[i+1]-Theta1[i-1])/(2*h):d8:=(Theta1[i+1]-2*Theta1[i]+Theta1[i-1])/h^2:

#Discretization Scheme
for i to N do

    eq1[i]:= (a1+6*delta1*(d1)^2)*(d2)-a2*R*(d1)+lambda0-a5*M^2*U0[i]+la*R*a2*Theta0[i]*(1+Qc*Theta0[i]) = 0:
    eq2[i]:= (a1+g*l+6*delta1*(d1)^2)*(d4)+12*delta1*(d1)*(d3)*(d2)-a2*R*(d3)-(H^2*a2*l+M^2*a5)*U1[i]+lambda1+la*R*a2*Theta1[i]*(1+2*Qc*Theta0[i])=0:
    eq3[i]:= (a4+(4*Rd*(1/3))*(1+(Tp-1)^3*Theta0[i]^3+3*(Tp-1)^2*Theta0[i]^2+(3*(Tp-1))*Theta0[i]))*(d6)/Pr-a3*R*(d5)+a5*Ec*M^2*U0[i]^2+a1*Ec*d1^2+2*Ec*delta1*d1^4+4*Rd*(Tp-1)*d5^2*(1+(2*(Tp-1))*Theta0[i]+(Tp-1)^2*Theta0[i]^2)/Pr= 0:
eq4[i]:= (a4+(4*Rd*(1/3))*(1+(Tp-1)^3*Theta0[i]^3+3*(Tp-1)^2*Theta0[i]^2+(3*(Tp-1))*Theta0[i]))*d8/Pr-a3*R*d7-l*a3*H^2*Theta1[i]+2*a5*Ec*M^2*U0[i]*U1[i]+2*a1*Ec*d1*d3+8*Ec*delta1*d1^3*d3+4*Rd*(Tp-1)*(Theta1[i]+(2*(Tp-1))*Theta0[i]*Theta1[i]+(Tp-1)^2*Theta0[i]^2*Theta1[i])*d6/Pr+8*Rd*(Tp-1)*d5^2*((Tp-1)^2*Theta0[i]*Theta1[i]+(Tp-1)*Theta1[i])/Pr+8*Rd*(Tp-1)*d5*d6*(1+(2*(Tp-1))*Theta0[i]+(Tp-1)^2*Theta0[i]^2)/Pr=0;
end do:

>	`union`( (seq( indets( eq1[i], name), i=1..N)), (seq( indets( eq2[i], name), i=1..N)), (seq( indets( eq3[i], name), i=1..N)), (seq( indets( eq4[i], name), i=1..N)) ); numelems(%):

{Theta0[1], Theta0[2], Theta0[3], Theta0[4], Theta0[5], Theta0[6], Theta0[7], Theta0[8], Theta0[11], Theta1[1], Theta1[2], Theta1[3], Theta1[4], Theta1[5], Theta1[6], Theta1[7], Theta1[8], Theta1[11], U0[1], U0[2], U0[3], U0[4], U0[5], U0[6], U0[7], U0[8], U0[11], U1[1], U1[2], U1[3], U1[4], U1[5], U1[6], U1[7], U1[8], U1[11]}

(3)

>

>	fsolve( eval( [ (seq(eq1[i], i=1..N)), (seq(eq2[i], i=1..N)), (seq(eq3[i], i=1..N)), (seq(eq4[i], i=1..N)) ] ) );

(4)

>

Download FDM_practice.mw

Implementation of finite difference and ...

Answered: KIRAN SAJJAN 40

May 04 2023

0 0

Please help me me to sosolve the equations with FDM i want to compair the both answer

i have already given the d-solve equations it is exicuting but FDM equations getting error i want answer from both the methods.

>

restart; with(plots)

>

a1 := 1/(1-p1-p2)^2.5

>

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

>

a4 := ((ks1*p1+ks2*p2)/(p1+p2)+2*kf+2*(ks1*p1+ks2*p2)-2*(p1+p2)*kf)/((ks1*p1+ks2*p2)/(p1+p2)+2*kf-ks1*p1-ks2*p2+(p1+p2)*kf)

>

a5 := 1+3*((p1*sigma1+p2*sigma2)/`σf`-p1-p2)/(2+(p1*sigma1+p2*sigma2)/((p1+p2)*`σf`)-((p1*sigma1+p2*sigma2)/`σf`-p1-p2))

>

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

>

interface(rtablesize = 100); interface(displayprecision = 5); Matrix([[Y, Nu1a, `Ν1b`, Nu1c], seq([j, seq([(a4+4*Rd*(1/3))*((eval(diff(Theta0(Y), Y), Ans[k]))(j)+0.1e-2*evalf(exp(l*t))*(eval(diff(Theta1(Y), Y), Ans[k]))(j))][], k = 1 .. numelems(RVals))], j = 0)]); interface(rtablesize = 10); interface(displayprecision = -1)

Matrix(%id = 18446746192273521350)

(1)

>

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

plotA:= display
  ( [ seq
      ( odeplot
        ( Ans[k],[Y,U0(Y)],
          Y=0..1,
          color=cols[k]
        ),
        k=1..numelems(RVals)
      )
    ],
    'axes'= 'boxed',labels=[Y,'U[0]'],title= Typesetting:-mtext("Steady Velocity distribution" ,color= Red)
  );
plotB:= display
  ( [ seq
      ( odeplot
        ( Ans[k],[Y,U1(Y)],
          Y=0..1,
          color=cols[k]
        ),
        k=1..numelems(RVals)
      )
    ],
    'axes'= 'boxed',labels=[Y,'U[1]'],title= Typesetting:-mtext("Unsteady Velocity distribution" ,color= Blue)

  );

>

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

plotE:= display( [ seq( odeplot
        ( Ans[k],[Y,Theta0(Y)],
          Y=0..1,
          color=cols[k]
        ),
        k=1..numelems(RVals)
      )
    ],
    'axes'= 'boxed',labels=[Y,'Theta[0]'],title= Typesetting:-mtext("Steady Temperature distribution" ,color= Red)
  );
plotF:= display
  ( [ seq
      ( odeplot
        ( Ans[k],[Y,Theta1(Y)],
          Y=0..1,
          color=cols[k]
        ),
        k=1..numelems(RVals)
      )
    ],
    'axes'= 'boxed',labels=[Y,'Theta[1]'],title= Typesetting:-mtext("Unsteady Temperature distribution" ,color= Blue)

  );

>

d1:= (U0[i+1]-U0[i-1])/(2*h):
d2:=(U0[i+1]-2*U0[i]+U0[i-1])/h^2:
d3:= (U1[i+1]-U1[i-1])/(2*h):d4:=(U1[i+1]-2*U1[i]+U1[i-1])/h^2:
d5:= (Theta0[i+1]-Theta0[i-1])/(2*h):d6:=(Theta0[i+1]-2*Theta0[i]+Theta0[i-1])/h^2:
d7:= (Theta1[i+1]-Theta1[i-1])/(2*h):d8:=(Theta1[i+1]-2*Theta1[i]+Theta1[i-1])/h^2:
h := 0.5:
ode1 := (a1+6*delta1*(d1)^2)*(d2)-a2*R*(d1)+lambda0-a5*M^2*U0[i]+la*R*a2*Theta0[i]*(1+Qc*Theta0[i]) = 0:
ode2 := (a1+g*l+6*delta1*(d1)^2)*(d4)+12*delta1*(d1)*(d3)*(d2)-a2*R*(d3)-(H^2*a2*l+M^2*a5)*U1[i]+lambda1+la*R*a2*Theta1[i]*(1+2*Qc*Theta0[i])=0:
ode3 := (a4+(4*Rd*(1/3))*(1+(Tp-1)^3*Theta0[i]^3+3*(Tp-1)^2*Theta0[i]^2+(3*(Tp-1))*Theta0[i]))*(d6)/Pr-a3*R*(d5)+a5*Ec*M^2*U0[i]^2+a1*Ec*d1^2+2*Ec*delta1*d1^4+4*Rd*(Tp-1)*d5^2*(1+(2*(Tp-1))*Theta0[i]+(Tp-1)^2*Theta0[i]^2)/Pr= 0:
ode4 := (a4+(4*Rd*(1/3))*(1+(Tp-1)^3*Theta0(Y)^3+3*(Tp-1)^2*Theta0(Y)^2+(3*(Tp-1))*Theta0(Y)))*(diff(Theta1(Y), Y, Y))/Pr-a3*R*(diff(Theta1(Y), Y))-l*a3*H^2*Theta1(Y)+2*a5*Ec*M^2*U0(Y)*U1(Y)+2*a1*Ec*(diff(U0(Y), Y))*(diff(U1(Y), Y))+8*Ec*delta1*(diff(U0(Y), Y))^3*(diff(U1(Y), Y))+4*Rd*(Tp-1)*(Theta1(Y)+(2*(Tp-1))*Theta0(Y)*Theta1(Y)+(Tp-1)^2*Theta0(Y)^2*Theta1(Y))*(diff(Theta0(Y), Y, Y))/Pr+8*Rd*(Tp-1)*(diff(Theta0(Y), Y))^2*((Tp-1)^2*Theta0(Y)*Theta1(Y)+(Tp-1)*Theta1(Y))/Pr+8*Rd*(Tp-1)*(diff(Theta0(Y), Y))*(diff(Theta1(Y), Y))*(1+(2*(Tp-1))*Theta0(Y)+(Tp-1)^2*Theta0(Y)^2)/Pr=0:

>

eta[0]:=0:
U0[0]:=0:U1[0]:=0:
Theta0[0] := 0:Theta1[0] := 0:
U0[1]:=0:U1[1]:=0:
Theta0[1] := 1:Theta1[1] := 0:
for i from 0 to 9 do
    eta[i+1]:=eta[i]+h;
    Eqa[i]:=eval(ode1,eta=eta[i]);
    Eqb[i]:=eval(ode2,eta=eta[i]);
    Eqc[i]:=eval(ode3,eta=eta[i]);
    Eqd[i]:=eval(ode4,eta=eta[i]);
end do:

>

Download FDM_practice.mw

Error in ANN...

Answered: KIRAN SAJJAN 40

April 28 2023

0 0

@Rouben Rostamian i want to implement for the below mentioned equations getting an error

>	restart; with(plots):

>	A := 0.5: M := 0.5: R := 0.5: Pr := 0.72: Ec := 0.1: g := 0.1:

>

ODE := diff(f(eta), eta, eta, eta) + f(eta)*diff(f(eta), eta, eta) - diff(f(eta), eta)^2 - M*diff(f(eta), eta) - A*diff(f(eta), eta) - A*eta*diff(f(eta), eta, eta)/2 = 0, (1 + R)*diff(Theta(eta), eta, eta) + Pr*(f(eta)*diff(Theta(eta), eta) - Theta(eta)*diff(f(eta), eta)) - Pr*A*(Theta(eta) + eta/2*diff(Theta(eta), eta)) + Pr*(Ec*diff(f(eta), eta, eta)^2 + g*Theta(eta)):
BCS:=f(0)=0,D(f)(0)=1,Theta(0)=1,D(f)(10)=0,D(Theta)(10)=0:
sol:=dsolve({ODE,BCS},numeric,output = listprocedure):
p1:=odeplot(sol,[eta,diff(f(eta),eta)],0..5,color=green,linestyle=3,thickness=3,legend = ["dsolve"]);
p2:=odeplot(sol, [eta, Theta(eta)], 0 .. 5, color = green, linestyle = 3, thickness = 3, legend = ["dsolve"]);

>	fVals:= Matrix([[eta, "f''(eta)"], seq( [j, eval(diff(f(eta),eta,eta), sol)(j)], j=0)]); ThetaVals:=Matrix([[eta, "Theta'(eta)"], seq( [j, eval(-diff(Theta(eta),eta), sol)(j)], j=0)]);

Matrix(2, 2, {(1, 1) = eta, (1, 2) = "f''(eta)", (2, 1) = 0, (2, 2) = -1.3662374845926935})

Matrix(%id = 18446746559113779070)

(1)

>

restart;

>

Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system

Error, (in plot) two lists or Vectors of numerical values expected

Download ANN.mw

FINITE ELEMENT METHOD...

Answered: KIRAN SAJJAN 40

April 21 2023

0 0

How to impliment finite element method of the below equations

restart;

OdeSys := diff(U(Y), Y, Y)+Theta(Y)+N*(Theta(Y)*Theta(Y))-(M*M)*U(Y) = 0, diff(Theta(Y), Y, Y)+E*(diff(U(Y), Y))^2 = 0;

Cond := U(0) = lambda*(D(U))(0), Theta(0) = A+g*(D(Theta))(0), U(1) = 0, Theta(1) = B; sys := [OdeSys, Cond];

Ans := dsolve(sys);

rectify my error sir...

Answered: KIRAN SAJJAN 40

March 08 2023

0 0

>

$restart; ContoursWithLabels := proc (Expr::algebraic, Range1::(name = range(realcons)), Range2::(name = range(realcons)), { contours::{posint, ({list, set})(realcons)} := 8, ImplicitplotOptions::(({list, set})({name, name = anything})) := NULL, GraphicOptions::(({list, set})({name, name = anything})) := NULL, TextOptions::(({list, set})({name, name = anything})) := NULL, Coloring::(({list, set})({name, name = anything})) := NULL }) local r1, r2, f, L1, h, S1, P, r, M, C, T, p, p1, m, n, i, x, y; x := lhs(Range1); y := lhs(Range2); f := unapply(Expr, x, y); if contours::posint then r1 := rand(convert(rhs(Range1), float)); r2 := rand(convert(rhs(Range2), float)); L1 := select(type, `~`[`@`(f, op)]({seq(([r1, r2])(), i = 1 .. 205)}), realcons); h := (L1[-6]-L1[1])/contours; S1 := [seq(L1[1]+(1/2)*h+h*(n-1), n = 1 .. contours)] else S1 := [contours[]] end if; userinfo(1, ContoursWithLabels, print('Contours' = evalf[2](S1)), `\n`); r := proc (k) options operator, arrow; rand(20 .. k-20) end proc; for C in S1 do P := plots:-implicitplot(Expr = C, Range1, Range2, gridrefine = 3, ImplicitplotOptions[]); for p in [plottools:-getdata(P)] do p1 := convert(p[3], listlist); n := nops(p1); if n < 500 then m := `if`(40 < n, (r(n))(), round((1/2)*n)); M[`if`(40 < n, [p1[1 .. m-11], p1[m+11 .. n]], [p1])[]] := NULL; T[[p1[m][], evalf[2](C)]] := NULL else h := trunc((1/2)*n); m := (r(h))(); M[p1[1 .. m-11], p1[m+11 .. m+h-11], p1[m+h+11 .. n]] := NULL; T[[p1[m][], evalf[2](C)], [p1[m+h][], evalf[2](C)]] := NULL end if end do end do; plots:-display([`if`(Coloring = NULL, NULL, plots:-densityplot(Expr, Range1, Range2, Coloring[])), plot([indices(M, 'nolist')], color = black, GraphicOptions[]), plots:-textplot([indices(T, 'nolist')], TextOptions[])], 'axes' = 'box', 'gridlines' = false, _rest) end proc$

>

fcns := { T(eta), f(eta)}; bc1 := a1=0.8578216437,a2=1.507519356,a3=1.191489352,a4=0.9767574138,a5=1.251405398,M=1, kp=1, Pr=6.8, N=1, g=0.5, A=1, B=0, Ec=1, lambda=0.5;
sys := (diff(f(eta), eta, eta))/(a1*a2)+T(eta)+N*(T(eta)*T(eta))-a3*(M*M)*f(eta)/a2-(kp*kp)*f(eta)/(a1*a2)= 0, a5*(diff(T(eta), eta, eta))/a4+Pr*Ec*((diff(f(eta), eta))^2+f(eta)^2*(kp*kp))/(a1*a2)= 0;
bc := f(0) = lambda*(D(f))(0), T(0) = A+g*(D(T))(0), f(1) = 0, T(1) = B;
R := dsolve(eval({bc1,bc, sys}), fcns, type = numeric, method = bvp[midrich], output = operator):
psi1 := unapply(eval(x*f(eta),R), x, eta);
psi2 := unapply(eval(x*T(eta),R), x, eta);

(diff(diff(f(eta), eta), eta))/(a1*a2)+T(eta)+N*T(eta)^2-a3*M^2*f(eta)/a2-kp^2*f(eta)/(a1*a2) = 0, a5*(diff(diff(T(eta), eta), eta))/a4+Pr*Ec*((diff(f(eta), eta))^2+f(eta)^2*kp^2)/(a1*a2) = 0

Error, (in dsolve/numeric/process_input) invalid specification of initial conditions, got {A = 1, B = 0, Ec = 1, M = 1, N = 1, Pr = 6.8, a1 = .8578216437, a2 = 1.507519356, a3 = 1.191489352, a4 = .9767574138, a5 = 1.251405398, g = .5, kp = 1, lambda = .5}

Error, invalid input: eval received R, which is not valid for its 2nd argument, eqns

>

ContoursWithLabels(
     psi1(eta, x), eta = 0 .. 1, x = -5 .. 5, contours=12,
     Coloring= [colorstyle= HUE, colorscheme= ["Red","Yellow","Green","Blue"], style= surface],
     TextOptions= [font= [HELVETICA,BOLD,7], color= black],
     ImplicitplotOptions= [gridrefine= 4],

     size= [600,600]
);

Error, (in ContoursWithLabels) invalid subscript selector

>

ContoursWithLabels(psi2(eta, x), eta = 0 .. 1, x = -5 .. 5, contours = 12, Coloring = [colorstyle = HUE, colorscheme = ["Red", "Yellow", "Green", "Blue"], style = surface], TextOptions = [font = [HELVETICA, BOLD, 7], color = black], ImplicitplotOptions = [gridrefine = 4], size = [600, 600])