Gabriel samaila

35 Reputation

4 Badges

6 years, 203 days

MaplePrimes Activity


These are replies submitted by Gabriel samaila

@Carl Love 

thank you so much for your response, is now running perfectely. However i notice that there is a problem of convergency. I think that is why the author ishak used shooting method before solving it. I dont know if you have idea of shooting method?

In his paper, the last figure is round (complete circle) and in the present code, if  i extent eta to any thing more than ten, the solution will not converge but in Ishark work, the solution converger. Or do you think of other ways? Please help because i want to extended his work.

new.mw

 
 

Download new.mw

 

@Carl Love 

Thank you for attending to my problem

attach is the code you requested for, possibly i did not implement it as instructed. please see the code for more details on what i want to do (written in upper case with red colour)

new.mw

 

@Carl Love 

Thank you for attending to my problem

attach is the code you requested for, possibly i did not implement it as instructed. What actually i am trying to do is to plot f'' against S(episolon) as in Ishak paper. i have already ploted figure 1 (theta against eta) but i could not plot f'' against episolon

ishak.mw          

Link to paper abstract
[article pdf uploaded without copyright permission, removed by moderator]

@Carl Love 

thank you so much for your response, you alway respond to my request i really thank you. However, i could not run it using maple 16, i may have to look for the recent version and try it.

@Carl Love 

Thank you for your response. However  I tried to do as you suggested but i couldnt get it right. 

Attach is your code you ever sent me to solve a particular problem. i have modify it to solve Reddy et al paper which i want to extend. Could you please show me how to get the numeric data of those plots?.

thank you in antipation 

reddy.mw
Download reddy.mw

@tomleslie 

Thanks for your response I will try and another way I think is to use print instead of plot in maple then I get the .dat then i can import it into MATLAB. 

Is that possible?

@tomleslie 

you always help me sort every sort of problem on maple, i'm very graful for all your post. i just wish i'm smart like you.

i really appreciate your effort.

@tomleslie 

thank you so much for your response is very helpful,  i think  the code is not supporting Maple 2016 version, see the error massage. Also is it possible to write boundary condition for u0,w0,f0 and theta0 different from ui,wi,fi and theta i for i=1..N?

thank you
 

  restart;
  PDEtools[declare](f(x),theta(x),u(x),w(x), prime=x):
  N:= 4;

f(x)*`will now be displayed as`*f

 

theta(x)*`will now be displayed as`*theta

 

u(x)*`will now be displayed as`*u

 

w(x)*`will now be displayed as`*w

 

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

 

4

(1)

  F:=      x -> local i;
                add
                ( p^i*f[i](x),
                  i=0..N
                );
  Theta:=  x -> local i;
                add
                ( p^i*theta[i](x),
                  i = 0..N
                );
  U:=      x -> local i;
                add
                ( p^i*u[i](x),
                  i = 0..N
                );
  W:=      x -> local i;
                add
                ( p^i*w[i](x),
                  i = 0 .. N
                );

Error, reserved word `local` unexpected

 

  HPMEq := collect
           ( (1-p)*diff(F(x), x$2)+p*(diff(F(x), x$2)-k1*diff(F(x), x)-k2*F(x)),
              p
           );
  HPMEr := collect
           ( (1-p)*diff(Theta(x),x$2)+p*(diff(Theta(x),x$2)-k11*diff(Theta(x),x)+k12*diff(U(x),x)^2+k13*diff(W(x),x)^2+k14*Theta(x)),
              p
           );
  HPMEs := collect
           ( (1-p)*diff(U(x),x$2)+p*(diff(U(x),x$2)-R*diff(U(x),x)-A-k8*W(x)-k7*U(x)+k5*Theta(x)+k6*F(x)),
              p
           );
  HPMEt := collect
           ( (1-p)*diff(W(x),x$2)+p*(diff(W(x),x$2)-R*diff(W(x),x)+k9*U(x)-k10*W(x)),
              p
           );

(-k1*(diff(F(x), x))-k2*F(x))*p+diff(diff(F(x), x), x)

 

(-k11*(diff(Theta(x), x))+k12*(diff(U(x), x))^2+k13*(diff(W(x), x))^2+k14*Theta(x))*p+diff(diff(Theta(x), x), x)

 

(-R*(diff(U(x), x))-A-k8*W(x)-k7*U(x)+k5*Theta(x)+k6*F(x))*p+diff(diff(U(x), x), x)

 

(-R*(diff(W(x), x))+k9*U(x)-k10*W(x))*p+diff(diff(W(x), x), x)

(2)

#
# renamed equ[1] to equ1
#
  for i from 0 to N do
      equ1[i] := coeff(HPMEq, p, i) = 0;
  end do;
#
# renamed equa[1] to equ2
#
  for i from 0 to N do
      equ2[i] := coeff(HPMEr, p, i) = 0;
  end do;
#
# renamed equat[1] to equ3
#
  for i from 0 to N do
      equ3[i] := coeff(HPMEs, p, i) = 0;
  end do;
#
# renamed equati[1] to equ4
#
  for i from 0 to N do
      equ4[i] := coeff(HPMEt, p, i) = 0;
  end do

diff(diff(F(x), x), x) = 0

 

-k1*(diff(F(x), x))-k2*F(x) = 0

 

0 = 0

 

0 = 0

 

0 = 0

 

diff(diff(Theta(x), x), x) = 0

 

-k11*(diff(Theta(x), x))+k12*(diff(U(x), x))^2+k13*(diff(W(x), x))^2+k14*Theta(x) = 0

 

0 = 0

 

0 = 0

 

0 = 0

 

diff(diff(U(x), x), x) = 0

 

-R*(diff(U(x), x))-A-k8*W(x)-k7*U(x)+k5*Theta(x)+k6*F(x) = 0

 

0 = 0

 

0 = 0

 

0 = 0

 

Warning, inserted missing semicolon at end of statement

 

diff(diff(W(x), x), x) = 0

 

-R*(diff(W(x), x))+k9*U(x)-k10*W(x) = 0

 

0 = 0

 

0 = 0

 

0 = 0

(3)

#
# Sequentially solve the system(s)
#
  for i from 0 to N do
      cons:= f[i](-1) = 1, f[i](1) = 1,
             theta[i](-1)=1, theta[i](1)=1,
             u[i](-1)=1, u[i](1)=1,
             w[i](-1)=1, w[i](1)=1;
      assign(dsolve( [ equ1[i], equ2[i], equ3[i], equ4[i], cons ]));
  od:

Error, (in dsolve) found the following equations not depending on the unknowns of the input system: {f[0](-1) = 1, f[0](1) = 1, u[0](-1) = 1, u[0](1) = 1, w[0](-1) = 1, w[0](1) = 1, theta[0](-1) = 1, theta[0](1) = 1}

 

#
# What happens next??
#
# OP can access any of the "individual" functions,
# for example theta[2](x) or w[3](x), just by entering
# these expressions, as in
#
  collect(theta[2](x), x);
  collect(w[3](x), x);
#
# Or one can get the "sum" terms for these functions
# just by entering F(x), Theta(x), U(x) or W(x). These
# expresssions are a bit lengthy
#
# For example, U(x) is shown below
#
  collect(Theta(x), x);

theta[2](x)

 

w[3](x)

 

Theta(x)

(4)

 


 

Download HOM2.mw

@tomleslie 

i said may be because I check it every seems to be same but I’m surprised is not coming up.

@tomleslie 

thank you so much for all your comment, they are all helpful, i have attach the source of the governing equation probalily i did not input it correctely. The table is not include in the pdf but i found the table in journal book compile by the author. i have tried to attach the journal book but is very heavy it couldnt attach.Sc.docx


 

  restart;
  kernelopts(version);

`Maple 16.00, IBM INTEL NT, Feb 23 2012, Build ID 731223`

(1)

#
# Define ODEs
#
  ODEs:= {diff(phi(r), r, r)+(1)/(r)*diff(phi(r),r)-Re*Sc*v(r)*(diff(phi(r), r))+Nt*(diff(theta(r), r, r)+(diff(theta(r), r))/r)/Nb = 0,
         diff(v(r), r, r)+(diff(v(r), r))/r-(Ha^2/(1-eta)^2+1/r^2)*v(r)-Re*v(r)*(diff(v(r), r)) = 0,
         diff(theta(r), r, r)+(1)/(r)*diff(theta(r),r)+Ec*Pr*(diff(v(r), r)-v(r)/r)^2-Pr*Re*v(r)*(diff(theta(r), r))+Nb*(diff(theta(r), r))*(diff(phi(r), r))+Nt*(diff(theta(r), r))^2 = 0};

{diff(diff(phi(r), r), r)+(diff(phi(r), r))/r-Re*Sc*v(r)*(diff(phi(r), r))+Nt*(diff(diff(theta(r), r), r)+(diff(theta(r), r))/r)/Nb = 0, diff(diff(v(r), r), r)+(diff(v(r), r))/r-(Ha^2/(1-eta)^2+1/r^2)*v(r)-Re*v(r)*(diff(v(r), r)) = 0, diff(diff(theta(r), r), r)+(diff(theta(r), r))/r+Ec*Pr*(diff(v(r), r)-v(r)/r)^2-Pr*Re*v(r)*(diff(theta(r), r))+Nb*(diff(theta(r), r))*(diff(phi(r), r))+Nt*(diff(theta(r), r))^2 = 0}

(2)

#
# Define boundary conditions: note that these restrict
# the solution to the range eta < r < 1.0 (or possibly
# 1.0 < r < eta, depending on whether eta>1 or eta<1)
#
# In the parameter list provided by the OP, eta is
# hardwired to 0.5, so valid solution will be in the range
#
#  0.5 <= r <= 1.0
#
  BCs := { phi(1) = 0, phi(eta) = 1, theta(1) = 0,
           theta(eta) = 1, v(1) = 0, v(eta) = 1
         }:

#
# Define values for Ha and Sc which are used
# to generarte output table
#
  HaVals:=[0, 2, 4, 6]:
  ScVals:=[0.5, 1.0, 2.0]:
#
# Solve the ODE system for all pairs of values in
# the lists HaVals, ScVals, and park all solutions
# in a Maple Table()
#
  for j from 1 by 1 to numelems(HaVals) do
      for k from 1 by 1 to numelems(ScVals) do
        #
        # Generate parameter list
        #
          pVals:= [ eta = .5,       Ha = HaVals[j],
                    Sc = ScVals[k], Nt = .1,
                    Nb = .1,        Re = 1,
                    Ec = 0.01,    Pr = 10
                  ];
        #
        # Produce ODE solution for this parameter list
        #
          sol[j,k]:=dsolve( `union`( eval(ODEs, pVals),
                                     eval(BCs, pVals)
                                   ),
                             numeric,
                             output=listprocedure
                           );
      od:
  od:

#
# generate plots for theta(y), phi(r), v(r) for
# all Ha and Sc values.
#
# Note that theta(r) and v(r) don't seem to vary
# as Sc changes value (but phi(r) does)
#
  pfunc:= z-> plot
              ( [ seq
                  ( seq
                    ( eval(z(r), sol[i,j])(r),
                      i=1..numelems(HaVals)
                    ),
                    j=1..numelems(ScVals)
                  ) ],
                  r=0.5..1,
                  title=cat( convert(z, string),
                             "(r) vs r"
                           ),
                  titlefont=[times, bold, 24],
                  legend=[ seq
                           ( seq
                             ( typeset
                               ( "Ha= ", HaVals[i],
                                 ", Sc= ", ScVals[j]
                               ),
                               i=1..numelems(HaVals)
                             ),
                             j=1..numelems(ScVals)
                           )
                         ],
                  legendstyle=[font=[times, bold, 16]]
              ):
  pfunc(theta);
  pfunc(phi);
  pfunc(v);

 

 

 

#
# Define function for the output table. This could
# be any expression involving  phi(y), theta(y), v(y)
# so just for illustration, let's use
#
#   phi(y)*theta(y)/v(r)5
#
  func:=phi(r):
#
# Define the value of 'y' at which the output is
# required. The boundary conditions require that
# 0.5 < r < 1.0. For illustration, let's use
#
#   r=0.75
#
  argr:=0.75:
#
# Print header lines for output table
#
  printf
  ( "%30s\n %8s%12.6f%12.6f%12.6f\n\n",
    "Sc","Ha", seq
               ( ScVals[j],
                 j=1..numelems(ScVals)
               )
  );

  seq
  ( printf
    ( " %8d%12.6f%12.6f%12.6f\n%",
      HaVals[k],
      seq
      ( eval
        ( func, sol[k,j]
        )(argr),
        j=1..numelems(ScVals)
      )
    ),
    k=1..numelems(HaVals)
  );

                            Sc
       Ha    0.500000    1.000000    2.000000

        0    0.146498    0.161720    0.192380
        2    0.207675    0.218275    0.239489
        4    0.284420    0.290105    0.301407
        6    0.325064    0.328319    0.334773

 

 


 

Download new_odeProb6.mw

@tomleslie 

eventhough is not specified at what value of phi(y) is given but i think is at phi(0). 

 

@acer 

sorry about the attachemen, i have been trying to get the attached table using same code but i couldnt get it.Sc.docx


 

restart;
kernelopts(version);

`Maple 16.00, IBM INTEL NT, Feb 23 2012, Build ID 731223`

(1)

ODEs:= {diff(phi(y), y, y)-Re*Sc*v(y)*(diff(phi(y), y))+Nt*(diff(theta(y), y, y)+(diff(theta(y), y))/y)/Nb = 0,
         diff(v(y), y, y)+(diff(v(y), y))/y-(Ha^2/(1-eta)^2+1/y^2)*v(y)-Re*v(y)*(diff(v(y), y)) = 0,
         diff(theta(y), y, y)+Ec*Pr*(diff(v(y), y)-v(y)/y)^2-Pr*Re*v(y)*(diff(theta(y), y))+Nb*(diff(theta(y), y))*(diff(phi(y), y))+Nt*(diff(theta(y), y))^2 = 0};

{diff(diff(phi(y), y), y)-Re*Sc*v(y)*(diff(phi(y), y))+Nt*(diff(diff(theta(y), y), y)+(diff(theta(y), y))/y)/Nb = 0, diff(diff(v(y), y), y)+(diff(v(y), y))/y-(Ha^2/(1-eta)^2+1/y^2)*v(y)-Re*v(y)*(diff(v(y), y)) = 0, diff(diff(theta(y), y), y)+Ec*Pr*(diff(v(y), y)-v(y)/y)^2-Pr*Re*v(y)*(diff(theta(y), y))+Nb*(diff(theta(y), y))*(diff(phi(y), y))+Nt*(diff(theta(y), y))^2 = 0}

(2)

BCs := {phi(1) = 0, phi(eta) = 1, theta(1) = 0, theta(eta) = 1, v(1) = 0, v(eta) = 1};

{phi(1) = 0, phi(eta) = 1, v(1) = 0, v(eta) = 1, theta(1) = 0, theta(eta) = 1}

(3)

HaVals:=[1, 2, 5, 10]:
for j from 1 by 1 to numelems(HaVals) do
    pVals := [eta = .5, Ha = HaVals[j], Sc = .8, Nt = .1, Nb = .1, Re = 2, Ec = 0.1e-1, Pr = 10];
    sol[j]:=dsolve( `union`(eval(ODEs, pVals), eval(BCs, pVals)),
                     numeric,
                     output=listprocedure
                  );
od:
displ_opt:= titlefont=[times, bold, 20],
            legend=[seq(typeset("Ha= ", HaVals[j]),j=1..4)],
            legendstyle=[font=[times, bold, 16]],
            color=[red, green, blue, black]:
plot( [seq( eval(theta(y), sol[j])(y),j=1..4)],
      y=0.5..1,
      title="theta(y) vs y",
      displ_opt
    );
plot( [seq( eval(phi(y), sol[j])(y),j=1..4)],
      y=0.5..1,
      title="phi(y) vs y",
      displ_opt
     );
plot( [seq( eval(v(y), sol[j])(y),j=1..4)],
      y=0.5..1,
      title="v(y) vs y",
      displ_opt
   );

 

 

 

 

 


 

Download odeProb4.mw

@acer 

yes I’m using maple 2016 and code is running now

@tomleslie 

i said the code has been running in maple 2016 but I had to remove the command I have shown you above before it work. 

Thank you

@tomleslie 

thank you for your recent post, i'm using maple 2016. the code can run  perfectely without the command

size=[800,600],

1 2 3 4 Page 1 of 4