Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 181 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

No...

tomleslie 13876
December 05 2018
0 0

@666 jvbasha 

The only way I can think opf doing this is to use the plots:-shadebetween() command which allows you to shade the area between two curves. There doesn't seem to be any way to combine this with the odeplot() command, so one has to produce the desired plots in a slightly different way (whihc isn't difficult.

See the attached - as always plots render better in a worksheet than they do on this site

> 

  restart;
  with(plots):
  fcns := {T(eta), f(eta)}:
  ep:= .1: M := 1: kp := .5: n := 1:
  ec:= .1: pr := 1: s := .1: N := 5:
  sys:= diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta$2)-diff(f(eta),eta)*diff(f(eta),eta)+ep*ep+(M+1/kp)*(ep-diff(f(eta),eta)) = 0,
        diff(T(eta),eta$2)+pr*(f(eta)*diff(T(eta),eta)-n*diff(f(eta),eta)*T(eta))+pr*(ec*diff(f(eta),eta$2)*diff(f(eta), eta$2)+ec*(M+1/kp)*diff(f(eta),eta)^2+s*T(eta)) = 0:
  bc:= f(0)=0, D(f)(0)=1, D(f)(N) = ep, T(0) = 1, T(N) = 0:
  R:= dsolve(eval({bc, sys}), numeric, output=listprocedure):
  psi := [-0.9e-1, -0.7e-1, -0.4e-1, -0.1e-1, 0, 0.1e-1, 0.4e-1, 0.7e-1, 0.9e-1]:
  for i to 9 do
      plt[i]:=odeplot( R, [eta, psi[i]/f(eta)], eta=0..5, axes=boxed);
  end do:
  display( [seq(plt[i], i=1..9)]);

 
> 

cols:=[magenta, maroon, red, orange, coral, yellow, green]:
plts:=[seq( plots:-shadebetween( psi[j+1]/eval( f(eta), R)(eta),psi[j]/eval( f(eta), R)(eta), eta=0..5, color=cols[j]), j=1..7)]:
plots:-display(plts);

 
> 

 

Download odeProb2.mw

Or...

tomleslie 13876
December 05 2018
0 0

@brian bovril 

All three versions shown on the wiki page ;-)

> 

  restart:
> 

  with(Student[Calculus1]):
  assume(r>0, h>0, h<r):
  interface(showassumed=0):
  SurfArea:= SurfaceOfRevolution(sqrt(r^2-(x-r)^2), x=0..h);
  Volume:= VolumeOfRevolution(sqrt(r^2-(x-r)^2), x=0..h);

2*h*Pi*r

  

Pi*r*h^2-(1/3)*Pi*h^3

 (1)
> 

#
# r, a, and h are related by the expression
#
  expr:=r^2=a^2+(r-h)^2:
#
# So isolate() 'r' from expr, and substitute
# the result in the Area & Volume expressions
# obtained previously
#
  ans:=isolate( expr, r):
  surf:=subs( ans, SurfArea);
  vol:=simplify(subs( ans, Volume));

Pi*(a^2+h^2)

  

(1/6)*Pi*h*(3*a^2+h^2)

 (2)
> 

#
# r, h, and theta are related by the expression
#
  expr:=r-h=r*cos(theta);
#
# So isolate() 'h' from expr, and substitute
# the result in the Area & Volume expressions
# obtained previously
#
  ans:=isolate( expr, h):
  surf:=collect(subs( ans, SurfArea),r);;
  vol:=simplify(subs( ans, Volume));

r-h = r*cos(theta)

  

2*(-cos(theta)+1)*Pi*r^2

  

(1/3)*Pi*r^3*(cos(theta)-1)^2*(2+cos(theta))

 (3)
> 

 

Download sphericalCap2.mw

Clarification needed...

tomleslie 13876
December 05 2018
0 0

Your ode system contains the independent variable 'eta'.

Your boundary condition contain the function eta(N), which occurs nowhere else in the worksheet.

What do you mean by eta(N)?

My last try...

tomleslie 13876
December 05 2018
0 0

@Cavemaann 

You state

the final solution we got for W(x) was not in terms of sin, it was W(x) = D4*cos(n*Pi*x/L), and since D4 cannot be zero the cos term has to be zero but integers of n gives us the cos expression to be equal to 1 not zero :(

You may have got the above, but I didn't. I got W(x) = D1*cosh(n*Pi*x/L). There is a difference between cos() and cosh() which seems to have escaped you

You also state

why cant "n" be anything else but an integer, its no fun

You supplied the boundary conditions for the problem c- not me. Satisfying these boundary conditions  requires that alpha*L=n*Pi, with n being an integer. If you want something different, you change the boundary conditions!

Then you want

Anyway, if i wanted to consider a concentrated mass on the beam, how would i have to manipulate the PDE (govering equation) to consider that into account? Any jump start on that :D

I would advise against doing this. It is obvious that you do not understand the simplest guided beam problem. Untile you do undestannd the simple prolem, you really shouldn't try a more complicated one!

Attached is the last worksheet I intend producing on this matter. It is an amalgamation of the two worksheets I posted previously with a few new comments. I have also added an animation which shows the beam vibration for the first five modes. In order to produce this animation, I just "guessed" some values for all of the parameters of the problem


 

``

> 

  restart;
#
# Define PDE
#
  pde:= E*M*diff(w(x,t),x$4)+A*rho*diff(w(x,t),t$2)=0;

E*M*(diff(diff(diff(diff(w(x, t), x), x), x), x))+A*rho*(diff(diff(w(x, t), t), t)) = 0

 (1)
> 

  test1:=w(x,t)=W(x)*cos(omega*t);
  pde1:=expand(eval(pde, test1)/(cos(omega*t)*E*M));
#
# Tidy up all the constants
#
  ode:=algsubs(A*rho*omega^2/(E*M)=alpha^4, pde1);
  sol:= dsolve(ode);
  sol:= convert(sol, trig);
############################################################
# Op seems(?) to have an issue because the simplify()
# command below does not "collect" terms in sinh() and
# cosh() in Maple 15 - so I can think of two alternatives
# which *might* work in Maple 15.
#
# Easiest way to do this is to define three possible
# simplification functions
#############################################################

  simp1:= z-> simplify(z):
  simp2:= z-> collect( z,
                       [ sinh(alpha*x),
                         cosh(alpha*x),
                         sin(alpha*x),
                         cos(alpha*x)
                       ]
                     ):
  simp3:= z->  collect
               ( collect
                 ( collect
                   ( collect
                     ( sol,
                       sinh(alpha*x)
                     ),
                     cosh(alpha*x)
                   ),
                   sin(alpha*x)
                 ),
                 cos(alpha*x)
               ):
#########################################################
# Now *one* of the following *ought* to work in Maple 15
# All three work in Maple 2018. Make sure you only
# uncomment one of the following - one of them should
# work!
########################################################
 
#  sol:=simp1(sol);
#  sol:=simp2(sol);
   sol:=simp3(sol);

  sol:=subs( [ op([2,1,1],sol)=D1,
               op([2,2,1],sol)=D2,
               op([2,3,1],sol)=D3,
               op([2,4,1],sol)=D4
             ],
             sol
           );
#
# Try the first couple of boundary conditions
#
  b1:= simplify(rhs(eval( diff(sol, x), x=0))=0);
  b2:= simplify(rhs(eval( diff(sol, x$3), x=0))=0);
#
# Modified the following slightly to avoid hard-
# wiring in the names of the solution variables
#
  bc1:= solve( [b1, b2],
               indets(b1, 'name') minus {alpha}
             );
  sol:= eval(sol, bc1);
 

w(x, t) = W(x)*cos(omega*t)

  

diff(diff(diff(diff(W(x), x), x), x), x)-A*rho*W(x)*omega^2/(E*M) = 0

  

-W(x)*alpha^4+diff(diff(diff(diff(W(x), x), x), x), x) = 0

  

W(x) = _C1*exp(-alpha*x)+_C2*exp(alpha*x)+_C3*sin(alpha*x)+_C4*cos(alpha*x)

  

W(x) = _C1*(cosh(alpha*x)-sinh(alpha*x))+_C2*(cosh(alpha*x)+sinh(alpha*x))+_C3*sin(alpha*x)+_C4*cos(alpha*x)

  

W(x) = (_C1+_C2)*cosh(alpha*x)+(-_C1+_C2)*sinh(alpha*x)+_C3*sin(alpha*x)+_C4*cos(alpha*x)

  

W(x) = D1*cosh(alpha*x)+D2*sinh(alpha*x)+D3*sin(alpha*x)+D4*cos(alpha*x)

  

alpha*(D2+D3) = 0

  

alpha^3*(D2-D3) = 0

  

{D2 = 0, D3 = 0}

  

W(x) = D1*cosh(alpha*x)+D4*cos(alpha*x)

 (2)
> 

          

#
# Now try the second pair of boundary conditions
#
  b3:= simplify(rhs(eval( diff(sol, x), x=L))/alpha=0);
  b4:= simplify(rhs(eval( diff(sol, x$3), x=L))/alpha^3=0);
#
# Consider
#
  b4-b3;
#
# This can only be true if either
#
#      D2=0, or
#      alpha*L=n*Pi (n an integer)
#
# Now consider
#
  b4+b3;
#
# This can only be true if either
#
#      D4=0, or
#      alpha*L=0
#
# I think(?) the only combination of the above two
# conditions which leads to non-trivial solutions
# is to have D4=0, and alpha*L=n*Pi. The latter
# condition however means restricting one of the
# variables in the definition of alpha - maybe
# omega??? Consider this later - for now, just make
# the substitutions D4=0, and alpha=n*Pi/L
#
  sol:=eval(sol, [D4=0, alpha=n*Pi/L]);
 

D1*sinh(alpha*L)-D4*sin(alpha*L) = 0

  

D1*sinh(alpha*L)+D4*sin(alpha*L) = 0

  

2*D4*sin(alpha*L) = 0

  

2*D1*sinh(alpha*L) = 0

  

W(x) = D1*cosh(n*Pi*x/L)

 (3)
> 

#
# Rebuild the complete solution
#
  ans:=eval(test1, sol);

w(x, t) = D1*cosh(n*Pi*x/L)*cos(omega*t)

 (4)
> 

################################################
#
# Get the natural frequencies of the beam
#
################################################
# Notice that earlier, alpha=n*Pi/L was
# required, which in turn means that there
# is a restriction on omega, which can be
# determined as
#
  omega__n:= rhs~
             ( [ ( allvalues
                   ( isolate
                     ( A*rho*omega^2/(E*M)=(n*Pi/L)^4,
                       omega
                     )
                   )
                 )
               ]
             );
#
# Assuming negative frequencies are not allowed
# substitute the first of the omega__n values
# (ie the positive one) in the general solution.
#
# In other words the modal frequencies omega__n
# are given by
#
#    omega__n = +sqrt(n^4*Pi^4*E*M/(L^4*A*rho))
#
# with n=1, 2, 3, 4, 5
#
# First output below shows the natural frequencies
# Second output shows the complete solution with
# the expression for the natural frequencies
# inserted appropriately
#
  ans:=eval(ans, omega=omega__n[1]);

[(n^4*Pi^4*E*M/(L^4*A*rho))^(1/2), -(n^4*Pi^4*E*M/(L^4*A*rho))^(1/2)]

  

w(x, t) = D1*cosh(n*Pi*x/L)*cos((n^4*Pi^4*E*M/(L^4*A*rho))^(1/2)*t)

 (5)
> 

#
# In order to plot the oscillations of this beam,
# insert some random set of values for all of the
# constants in the problem
#
  pltexpr:= rhs
            ( eval
              ( algsubs
                ( E*M/(A*rho)=1, ans ),
                [D1=1e-03, L=1]
              )
            );
#
# Produce an animation showing the first five
# vibration modes of the beam. Note that from
# the form of the answer obtained, one would
# expect the following
#
# 1, As n increases, the argument of the cos()
#    term increases: in other words the oscillation
#    frequency increases. So n=1 is a "slow"
#    vibration and n=5 is a "fast" vibration
#
# 2. As n increases, the argument of the cosh()
#    term increases: in other words the amplitude
#    of the oscillation increases (really rapidly!)
#
# Because of (2) above, the amplitudes of the
# various modes is arbitrarily scaled in the
# following animation, just so that everything
# can be seen on the same graph
#
#
  plots:-animate( plot,
                  [ [ eval(pltexpr, n=1),
                      eval(pltexpr, n=2)/20,
                      eval(pltexpr, n=3)/600,
                      eval(pltexpr, n=4)/12000,
                      eval(pltexpr, n=5)/240000
                    ],
                    x=0..1,
                    color=[red, green,blue, cyan, black],
                    legend=[ n=1, n=2, n=3, n=4, n=5],
                    legendstyle=[ font=[ times, bold, 20]]
                  ],
                  t=0..0.5,
                  frames=100
                  
                );

0.1e-2*cosh(n*Pi*x)*cos((n^4*Pi^4)^(1/2)*t)

  

 
> 

 

``


 

Download Beam2.mw

 

I should also add...

tomleslie 13876
December 04 2018
0 0

@tomleslie 

that performing the same calculation as before with the Grid:-Seq() command on a four-core machine, gives  ~3.5X speed-up: better than anything I saw in my earlier post using the Threads() package.

See the attached.

> 

  restart;
  kernelopts(version);
  kernelopts(multithreaded);
  numCP:=kernelopts(numcpus);
  N:=10^9:
  f:= proc(a,b) local j;
           add(j, j=a..b)
      end proc:
  CodeTools:-Usage(add([seq( f(k*N/4+1,(k+1)*N/4), k=0..3)]));
  CodeTools:-Usage(add([Grid:-Seq(f(k*N/4+1,(k+1)*N/4), k=0..3)]));

`Maple 2018.2, X86 64 WINDOWS, Nov 16 2018, Build ID 1362973`

  

true

  

8

  

memory used=2.50GiB, alloc change=-3.01MiB, cpu time=53.94s, real time=53.94s, gc time=2.54s

  

500000000500000000

  

memory used=1.16MiB, alloc change=19.69MiB, cpu time=78.00ms, real time=16.02s, gc time=31.20ms

  

500000000500000000

 (1)
> 

 

Download paraGrid.mw

Good luxk...

tomleslie 13876
December 04 2018
0 0

@acer 

I have been trying to explain to this OP for a while, that sin(alpha*L)=0 is satisfied whenever alpha*L=n*Pi, with n an integer.

However I now have to conclude that the OP does not understand the sin() function - pretty much puts him/her beyond help on this forum

I was trying to cover several points...

tomleslie 13876
December 04 2018
0 0

@acer 

distributed across a couple of posts here, namely

  1. Use of Threads:-Seq() does not necessarily improve real execution times -although playing with the 'tasksize' option can help
  2. There are two different methods (supplied by Maplesoft at different places in the help) for doing the same example calculation. One of these provides an improvement (2X), and one of these provides a substantial degradation (5X), and in the latter case I have absolutely no idea why!!

Basically a difference between Maple 15 ...

tomleslie 13876
December 04 2018
0 0

@Cavemaann 

because in my original worksheet I used a simplify() command which "combined" terms in sinh() (and cosh()). You state that this "combination" does not happen in Maple 15. I have to believe you: I can't check because I no longer have Maple 15 installed. Since I don't have access to Maple 15 any more, I cannot absolutely guarantee what will/will not work!

In the attached worksheet I have incorporated three different 'simplification' functions, named simp`1, simp2, simp3. Hopefully one of these will work correctly in Maple 15

  1. simp1() just a different way of applying the simplify() function, as used in my earlier worksheet - so this is probably not going to work for you!
  2. simp2() and simp3() are (slightly) different simplifications using the collect() command, which should 'force' collection of terms in sin(alpha*x), cos(alpha*x), sinh(alpha(x) and cosh(alpha*x)

You have to make sure to uncomment only one of these - make sure the other two function applications are commented out. Note all three of these simplification functions *work* in Maple 2018

In the attached, the application of simp1() and simp2() are commented out. I have highlighted code change reasons in red

For reasons I don't understand the worksheet will not display here, so you will have to download  to examine/execute/edit

Download Beam.mw

 

I'll try again...

tomleslie 13876
December 03 2018
0 0

@Cavemaann 

Don't provide pictures of code, provide code

In other words don't use the big green up-arrow to insert something un-executable and un-editable.

Use it to upload a complete executable/editable worksheet which demonstrates your problem

 

So much imprecision...

tomleslie 13876
December 03 2018
0 0

@Cavemaann 

You state that

we need to subsitute alpha*L = 4.73, 7.853, 10.996 to find the roots of the solution

quantiies on the rhs of the above expression cporrespond to (3/2)*Pi, (5/2)*Pi, (7/2)*Pi. In my original code I stated that

alpha*L=n*Pi, and thus alpha=n*Pi/L, with n any positive or negative integer

Please note that (3/2), (5/2), (7/2) are not integers

##########################################################################################

Impossible for me to diagnose the execution errors you are receiving if I cannot see the code you are running. Please run exactly the worksheet I provided in your Maple version and then post the resulting worksheet here using the big green up-arrow in the Mapleprimes toolbar.

If (as I suspect) my original code runs without error in your Maple version, then please post the worksheet you are actually running - the big green up-arrow is your friend

 

See the help page...

tomleslie 13876
December 02 2018
0 0

@mehdibaghaee 

at

?examples,Threads

Addendum...

tomleslie 13876
December 02 2018
0 0

I overlooked the fact that you were interested in the period of the oscillation. It can be a little awkward to calculate this efficiently from a numerical solution, so I have added an execution group in the attached to do this.

I also included a simple animation of the pendulum


 

> 

#
# Initialise
#
  restart:
#
# Define parameters and ODE
#
  L:=1: g:=9.81: omega:=sqrt(g/L): alpha:=0.25:
  de:=diff(theta(t),t$2)+2*alpha*diff(theta(t),t)+omega^2*sin(theta(t))=0;
#
# Make up some plausible(?) boundary conditions
# because OP doesn't provide any. So just set up
# for theta=Pi/2 (so pendulum is at the horizontal)
# with D(theta)=0 (so it is stationary). Then it
# is set free
#
  bcs:=theta(0)=Pi/2, D(theta)(0)=0:
#
# Solve the ODE numerically
#
  sol:=dsolve( {de,bcs},
               numeric,
               output=listprocedure
             ):
#
# plot the solution
#
  plots:-odeplot(sol, [t, theta(t)], t=0..10);

diff(diff(theta(t), t), t)+.50*(diff(theta(t), t))+9.810000002*sin(theta(t)) = 0

  

 
> 

#
# Produce a simple animation of the pendulum
# from the above solution
#
  with(plots):
  with(plottools):
  f:=eval(theta(t), sol):
  p3:=[ seq
        ( display
          ( [ line
              ( [0,0],
                [L*sin(f(t)), -L*cos(f(t))],
                thickness=5
              ),
              disk
              ( [L*sin(f(t)), -L*cos(f(t))],
                0.1,
                color=red
              )
            ]
          ),
          t=0..10, 0.1
        )
      ]:
  display( p3,
           axes=none,
           scaling=constrained,
           insequence=true
         );

 
> 

##########################################
# Determine period from numerical data.
##########################################
#
# Set up number of zero crossings to retain
# and initialise a matrix for the results.
# Columns in the matrix correspond to
#
#  col[1] = time value of zero crossing
#  col[2] = theta(t) value at "zero"-crossing.
#           Should be identically zero, but
#           in a numerical "rootfinding" process
#           some residual is inevitable
#  col[3] = slope of theta(t): will alternate
#           between '+' and '-' values,
#           corresponding to positive-going
#           and negative-going zero crosses
#  col[4] = period, based on subtracting
#           alternate col[1] values          
#
  nc:= 20:
  zc:= Matrix(nc+1, 4):
  zc[1, 1..4]:= < t, theta(t), diff(theta(t),t), period>:
#
# Set up dsolve() to output data whenever
# theta(t) crosses zero
#
  sol2:=dsolve( {de, bcs},
                numeric,
                events=[[theta(t),halt]]
              ):
#
# Suppress warnings generated by the 'events'
# process, and reset the size of the matrix
# to be displayed to accommodate the number
# of zero crossings requested
#
  interface(warnlevel=0, rtablesize=nc+1):
#
# Compute/store the zero-crossings
#
  for j from 2 by 1 to nc+1 do
      zc[j,1..3]:=<rhs~(sol2(50))[]>:
      if   j > 3
      then zc[j, 4]:= zc[j,1]-zc[j-2,1];
      fi:
      sol2(eventclear):
  od:
#
# Display output
#
  zc();
#
# Reset the warninglevel and rtablesize
# to defaults
#
  interface(warnlevel=3, rtablesize=10):

Matrix(21, 4, {(1, 1) = t, (1, 2) = theta(t), (1, 3) = diff(theta(t), t), (1, 4) = period, (2, 1) = .6269060985956958, (2, 2) = -0.5692061406e-16, (2, 3) = -3.839864958768343, (2, 4) = 0, (3, 1) = 1.7224713604826027, (3, 2) = 0.2426444462e-15, (3, 3) = 2.9548098407717602, (3, 4) = 0, (4, 1) = 2.7775834532815113, (4, 2) = 0.2882351476e-15, (4, 3) = -2.284128377963782, (4, 4) = 2.150677354685816, (5, 1) = 3.811821024282486, (5, 2) = -0.3583017130e-15, (5, 3) = 1.7700311601762109, (5, 4) = 2.089349663799883, (6, 1) = 4.834478807120099, (6, 2) = 0.3080353897e-15, (6, 3) = -1.373565368256029, (6, 4) = 2.0568953538385872, (7, 1) = 5.850470696990112, (7, 2) = 0.4160354786e-15, (7, 3) = 1.0667744542782682, (7, 4) = 2.0386496727076264, (8, 1) = 6.862546280552013, (8, 2) = 0.1402212222e-15, (8, 3) = -.8289074164082234, (8, 4) = 2.028067473431914, (9, 1) = 7.872293902560268, (9, 2) = 0.4898222127e-16, (9, 3) = .6442656300699325, (9, 4) = 2.021823205570156, (10, 1) = 8.880648201817461, (10, 2) = -0.4217563392e-15, (10, 3) = -.50084027610827, (10, 4) = 2.0181019212654485, (11, 1) = 9.888165176226751, (11, 2) = -0.1652324111e-15, (11, 3) = .3893847028506748, (11, 4) = 2.015871273666483, (12, 1) = 10.895177723427093, (12, 2) = -0.1949802082e-15, (12, 3) = -.3027512234814019, (12, 4) = 2.0145295216096315, (13, 1) = 11.90188594953875, (13, 2) = 0.8105088866e-16, (13, 3) = .2354016268713377, (13, 4) = 2.0137207733119986, (14, 1) = 12.90841040550298, (14, 2) = -0.1358098746e-15, (14, 3) = -.18303875439701034, (14, 4) = 2.0132326820758877, (15, 1) = 13.914823836798508, (15, 2) = 0.1122081486e-15, (15, 3) = .14232549767633665, (15, 4) = 2.012937887259758, (16, 1) = 14.92117015401817, (16, 2) = 0.6457440377e-16, (16, 3) = -.11066903342049417, (16, 4) = 2.0127597485151902, (17, 1) = 15.927475901756589, (17, 2) = -0.6501528442e-16, (17, 3) = 0.8605414634974683e-1, (17, 4) = 2.012652064958081, (18, 1) = 16.933757098767344, (18, 2) = 0.2645453301e-16, (18, 3) = -0.669143072256107e-1, (18, 4) = 2.0125869447491738, (19, 1) = 17.94002345162486, (19, 2) = -0.3098870038e-16, (19, 3) = 0.52031589797577334e-1, (19, 4) = 2.0125475498682714, (20, 1) = 18.946280819116208, (20, 2) = -0.4170959639e-16, (20, 3) = -0.4045906379533677e-1, (20, 4) = 2.012523720348863, (21, 1) = 19.952532707900076, (21, 2) = -0.2692633036e-16, (21, 3) = 0.314604689168979e-1, (21, 4) = 2.012509256275216})

 (1)
> 

 


 

Download pendulum2.mw

Check the graphs in the attached...

tomleslie 13876
December 01 2018
0 0

@Bagh Ali 

Note that these were obtained with a much hihger value of 'Pr' (ie 7.0) rather than the value you supplied (0.773). You can adjust this value just by changing the entry at

PrVal:=7.0:

whhic I have highlighted in red in the attached


 

> 

  restart;
  with(plots):
#
# Define the ODE system
#
  odeSys:= { diff(F(eta), eta$2)+2*eta*diff(F(eta),eta)-4*M^2*F(eta)=0,
             diff(theta(eta), eta$2)+2*Pr*eta*diff(theta(eta),eta)+Pr*Nb*diff(theta(eta), eta)*diff(phi(eta), eta)+Pr*Nt*diff(theta(eta),  eta)^2=0,
             diff(phi(eta), eta$2)+2*Pr*Le*eta*diff(phi(eta),eta)+(Nt/Nb)*diff(theta(eta), eta$2)=0
           }:
#
# Define the first set of boundary conditions
#
  bcs1:= { F(0)=1, theta(0)=1, phi(0)=1,
           F(inf)=0, theta(inf)=0, phi(inf)=0
         }:
#
# Need some parameter values for numeric
# solution but note that "infinity"=10
#
# Note also with Pr=0.773 (OP supplied value)
# graphs don't agree too well with those supplied
# by OP. Get better agreement Pr~7. Changing the
# setting below will allow this to be changed
# easily
#
  PrVal:=7.0:
  NtVals:=[0.07, 0.14, 0.21, 0.28]:
  for k from 1 by 1 to numelems(NtVals) do
      pList:=[ Le=0.5, M=0.5, Nb=0.1,
               Nt=NtVals[k], Pr=PrVal, inf=10
             ]:
      sol1[k]:= dsolve( eval
                        ( `union`( odeSys, bcs1),
                           pList
                       ),
                       numeric
                     );
     od:
  display
  ( [ seq
      ( odeplot
        ( sol1[i],
          [eta, theta(eta)],
          eta=0..2
        ),
        i=1..numelems(NtVals)
      )
    ],
    color = [red, green, blue, black],
    title = typeset( theta(eta), " versus ", eta),
    titlefont = [times, bold, 20]
  );

 
> 

  LeVals:=[0.1, 0.3, 0.5, 0.7]:
  for k from 1 by 1 to numelems(LeVals) do
      pList:=[ Le=LeVals[k], M=0.5, Nb=0.1,
               Nt=0.1,Pr=PrVal, inf=10
             ]:
      sol1[k]:= dsolve( eval
                        ( `union`( odeSys, bcs1),
                           pList
                       ),
                       numeric
                     );
     od:
  display
  ( [ seq
      ( odeplot
        ( sol1[i],
          [eta, phi(eta)],
          eta=0..2
        ),
        i=1..numelems(LeVals)
      )
    ],
    color = [red, green, blue, black],
    title = typeset( phi(eta), " versus ", eta),
    titlefont = [times, bold, 20]
  );

 
> 

 

NULL


 

Download bvivShoot2.mw

What ois Pr...

tomleslie 13876
December 01 2018
0 0

@Bagh Ali 

Your supplied ODE system contains the parameter 'Pr'. Nowhere in the additional data you have supplied is there a value for this parameter. Hence I cannot attempt to reproduce plots.

Well...

tomleslie 13876
December 01 2018
0 0

@Cavemaann 

  1. I got a "non-trivial" solution
  2. If you look at the execution group dealing with the second pair of boundary conditions, fulfilling the latter requires that D2*sin(alpha*L)=0 and D4*sinh(alpha*L)=0. There are various ways to fulfil these requirements.
    1. Dealing with the second requirement: sinh(x)=0 only if x=0 - so one possibility is that alpha*L=0: but this requires alpha=0, which results in trivial solutions everywhere. The other alternative is that D4=0
    2. in the first condition sin(x) will return zero whenever the argument x is an integer multiple of Pi, so alpha*L=n*Pi, and thus alpha=n*Pi/L, with n any positive or negative integer
  3. I can't produce any meaningful graphs unless the various constants (E, M, L, A, rho, phi) have numerical values. If you assign values to these quantiites then you can plot the final expression I obtained (ignoring the scaling factor D2) in a 3-D plot, as a function of 'x' and 't' parameterized by integer values of the quantity 'n'
First 79 80 81 82 83 84 85 Last Page 81 of 207