Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

Hello,

     Can you enter debug in a worksheet vs a proc?   

    I have a for loop that is crashing in the middle with a bad calculated number to provide to a thermoprops function.   Is there a way to single step through a for loop (1500 iterations) to see which line is getting tripped up?  I can not post the worksheet as it's private development.   The question is generic to any "for loop".

   I have also tried "porting" into a code edit region, but still can not enter the debugger mode.   Is proper debugger only available if I turn the "for loop" into a proc?   If yes, what is the availability of packages called right after restart for the worksheet for a proc, or do all the pre-calcs developed in the worksheet and the packages to be variables fed into the proc or packages re-declared insde a new proc ???

  Sorry if this is a trival question.

Thanks in advance,
Bill
 

I solve numerically a DAE system whose independent variable is named t and the dependent variables are d[1](t), ..., d[n](t).
I would like to to 2D or 3D plots of the solutions and color the resulting curve using a function f(...) of the remaining dependent variables.

Here is a simple example.

 

restart

with(plots):

0

(1)

sys := {
   diff(x(t), t) = v(t)
  ,diff(v(t), t) = cos(t)
  ,x(0) = -1
  ,v(0) = 0

  ,px(t) = piecewise(x(t) >=0, 1, -1)
  ,pv(t) = piecewise(v(t) >=0, 1, -1)
}:

sol := dsolve(sys, numeric):

 

odeplot(sol, [t, x(t), v(t)], t=0..4*Pi)
 

 

# I would like to color this space curve depending on the signs of x(t) and y(t)
#
# for instance, f being a "color function"
f := proc(s)
  local a, b:
  if s::numeric then
    a := round(eval(px(t), sol(s))):
    b := round(eval(pv(t), sol(s))):
    return piecewise(a+b = 2, "Green", a = 1, "Red", b = 1, "Blue", "Gold")
  end if:
end proc:

SOL := proc(s)
  if s::numeric then
    eval([t, x(t), v(t)], sol(s))
  end if:
end proc:


# I would like to make something like this to work

plot3d(SOL(s), s=0..4*Pi, colorfunc=f(s)):  #... which generates a void plot


 

# In some sense a continuous version of this

opts := symbol=solidbox, symbolsize=20:
display( seq(pointplot3d({SOL(s)}, opts, color=f(s)), s in [seq](0..6, 0.1)) );

 

 


 

Download coloring.mw


How can I fix (if possible) the syntax in the command 

plot3d(SOL(s), s=0..4*Pi, colorfunc=f(s)):

???

Thanks in advance

 

 

Has someone tried to connect Grasshopper with Maple yet?

Grasshopper is a visual programming environment on Rhino.

https://en.wikipedia.org/wiki/Grasshopper_3D

Hello;

I am facing "Error, too many levels of recursion" in loops. kindly guide me. Thanks

 

1D1P.mw


 

restart; printlevel := 3

restart; with(LinearAlgebra); with(linalg); with(CodeGeneration); with(plots); Digits := 30; `ε` := 0.1e-1; k := .5; m := 1; c := 1; q := 1

L[time] := 1:

.200000000000000000000000000000

 

.200000000000000000000000000000

 

.200000000000000000000000000000

(1)

``

ICff := 1.0*exp(-p^2/(2.0))*(1.0+epsilon*cos(1.0*k*x))/sqrt(2.0*Pi):

for i from 0 while i <= N[x] do for j from 0 while j <= N[p] do f[0, i, j] := eval(ICff, [x = i*`&Delta;x`, p = j*`&Delta;p`]) end do end do:

for n from 0 while n <= T do for j from 0 while j <= N[p] do f[n, 0, j] := 0; f[n, N[x], j] := 0 end do end do:

for n from 0 while n <= T do for i from 0 while i <= N[x] do f[n, i, 0] := 0; f[n, i, N[p]] := 0 end do end do:

``

for n1 from 0 while n1 <= T-1 do for i1 while i1 <= N[x] do for j1 while j1 <= N[p] do f[n1+1, i1, j1] := f[n1, i1, j1]-`&Delta;t`*j1*`&Delta;p`*(f[n1+1, i1+1, j1]-f[n1+1, i1-1, j1]+f[n1, i1+1, j1]-f[n1, i1-1, j1])/((4*`&Delta;x`)*(1.0)) end do end do end do

Error, too many levels of recursion

 


 

Maple 2021.2 gives error when calling int() first time. Second time it returns the integral unevaluated.

Is this a known issue? But it seems possible to trap this error for now,, which is good.
 

interface(version);

`Standard Worksheet Interface, Maple 2021.2, Windows 10, November 23 2021 Build ID 1576349`

Physics:-Version()

`The "Physics Updates" version in the MapleCloud is 1133 and is the same as the version installed in this computer, created 2022, January 20, 23:5 hours Pacific Time.`

restart;

w:=(7*x - 3 + sqrt(x^2 + (x^3*(x - 1)^2)^(1/3) - x) + sqrt(-2*((-x^2 + x + (x^3*(x - 1)^2)^(1/3)/2)*sqrt(x^2 + (x^3*(x - 1)^2)^(1/3) - x) + x^2*(x - 1))/sqrt(x^2 + (x^3*(x - 1)^2)^(1/3) - x)))/(12*x*(x - 1));

(1/12)*(7*x-3+(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2)+(-2*((-x^2+x+(1/2)*(x^3*(x-1)^2)^(1/3))*(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2)+x^2*(x-1))/(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2))^(1/2))/(x*(x-1))

int(w,x)

Error, (in IntegrationTools:-Indefinite:-AlgebraicFunction) invalid argument for sign, lcoeff or tcoeff

int(w,x)

int((1/12)*(7*x-3+(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2)+(-2*((-x^2+x+(1/2)*(x^3*(x-1)^2)^(1/3))*(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2)+x^2*(x-1))/(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2))^(1/2))/(x*(x-1)), x)

restart;
w:=(7*x - 3 + sqrt(x^2 + (x^3*(x - 1)^2)^(1/3) - x) + sqrt(-2*((-x^2 + x + (x^3*(x - 1)^2)^(1/3)/2)*sqrt(x^2 + (x^3*(x - 1)^2)^(1/3) - x) + x^2*(x - 1))/sqrt(x^2 + (x^3*(x - 1)^2)^(1/3) - x)))/(12*x*(x - 1));
try
  int(w,x);
catch:
  print("error happened");
end try;

(1/12)*(7*x-3+(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2)+(-2*((-x^2+x+(1/2)*(x^3*(x-1)^2)^(1/3))*(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2)+x^2*(x-1))/(x^2+(x^3*(x-1)^2)^(1/3)-x)^(1/2))^(1/2))/(x*(x-1))

"error happened"

 


 

Download int_problem_jan_21_2022.mw

Dear esteem Colleagues,

Please how do I modify the following two files (though similar) to get consistent errors? I am not sure where I made the mistake.

Any modifications would be appreciated.

Thank you all for your time and mentorship. Best regard

Biratu_Mapleprimes.mw

DDE_2_Mapleprime.mw

I would like to emphasize a particular region of my plot by shading the whole strip between x=0 and x=1. I learned of the command shadebetween, but the functions it takes as arguments for the boundaries of the region to shade must be functions of x.

Does someone know how to shade between two vertical lines?

I am a new user to Maple, but have used PTC's Mathcad and Prime.

I have several matrices that I need to combine both vertically and horizontally.  The command augment will combine the matrices horizontally.  I was unable to find a command to combine the matrices vertically.  I tried the stack command (Prime command).  The Maple stack command appears to have a different functionality.  Could you suggest which command would do the job?  The matrices are large and I would like to limit the number of rows and columns shown in the document.

The next step will be sorting the data. I am researching this at this time.

Thank you for the support,

David Tietje

Alternating serie

 

sum((-1)^(n+1)/(2*n-1), n = 1 .. infinity)

(1/4)*Pi

(1)

 

sum((-1)^(n + 1)/(2*n - 1), n = 1 .. infinity):

sum((-1)^(n+1)/(2*n-1), n = 1 .. 4)

76/105

(2)

expand(sum((-1)^(n + 1)/(2*n - 1), n = 1 .. 4),symbolic);

76/105

(3)

?expand

series(sum((-1)^(n + 1)/(2*n - 1), n = 1 .. infinity),n=0,5);

series((1/4)*Pi,n)

(4)

 

Info series

   

How to get for n= 4  "for (&sum;)(((-1)^(n+1))/(2 n-1))  =   symbolic term 1+ symbolic term 2+... "

 

 

Download onderzoek_reeks_-hoe_krijg_ik_een_partieke_symbolische_som.mw

Input:

 a := x^2;
 whattype(x);
 b := x[1]^2;
 whattype(x[1]);
 CodeGeneration[C](a);
 CodeGeneration[C](b);

Output:

Do you know why cg0 =/= x[0]*x[0]?

Hello people in MaplePrimes,

I can't understand a code with PLOT3D using GRID in it.
I hope you will give me explanation about it.

restart;
k:=[[1,2,10],[1,1,8],[0,1,5],[0,0,6]];
PLOT3D(GRID(1..2,1..3,k),
       AXESLABELS(x,y,z));

I cannot understand the correspondence between points in the graph and lists in the code such as [1,2,10], [1,1,8].

 As two ranges in GRID augments in the code are 1..2 and 1..3, x takes value from 1 to 2, and
y from 1 to 3 in the picture.
And, as k is k:=[[1,2,10],[1,1,8],[0,1,5],[0,0,6]], which means this list of lists are 4 X 3, so in the picture
as for x the interval of 1 to 2 is divided to three sections with the points of value of 0, 1/3, 2/3 and 3/3, which total number is 4.
And, as for y the interval of 1 to 3 is devided to two sections with the points of value of 1, 2 and 3, which total number is 3.

And, I cannot understand why from the grid of [1,1,8] in the code, which I think is the coodinate of (1, 1) on x-y plane, does not mean that the value of z on that point on the graph is not 8, but 1.

graph.mw

When I added kernelopts('assertlevel'=2): now Maple gives an error from call to solve. Also PDEtools:-Solve  gives same error. 

If I remove the  kernelopts('assertlevel'=2): it works OK and gives solution.

Is this known issue? I'd like to keep kernelopts('assertlevel'=2): and still use solve. Maple 2021.2 on windows 10

interface(version);

`Standard Worksheet Interface, Maple 2021.2, Windows 10, November 23 2021 Build ID 1576349`

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1131 and is the same as the version installed in this computer, created 2022, January 7, 9:5 hours Pacific Time.`

restart;

3

eq:=1/8*2^(1/2)*(-2*x*(w^4*x^11+(4*w^4-6*w^3)*x^10+(10*w^4-22*w^3+17*w^2)*x^9+(16*w^4-50*w^3+62*w^2-30*w)*x^8+(71/2+19*w^4-72*w^3+132*w^2-114*w)*x^7+(140-685/3*w+16*w^4-76*w^3+176*w^2)*x^6+(263-841/3*w+10*w^4-56*w^3+331/2*w^2)*x^5+(292-1385/6*w+105*w^2+4*w^4-30*w^3)*x^4+(821/4+w^4-733/6*w-10*w^3+91/2*w^2)*x^3+(173/2-2*w^3-122/3*w+11*w^2)*x^2+(77/4+3/2*w^2-37/6*w)*x+3/2-1/2*w)*(-1+I*3^(1/2))^(1/2)+2^(1/2)*(I*x*(w^4*x^11+(4*w^4-6*w^3)*x^10+(10*w^4-22*w^3+17*w^2)*x^9+(16*w^4-50*w^3+62*w^2-30*w)*x^8+(71/2+19*w^4-72*w^3+132*w^2-114*w)*x^7+(140-685/3*w+16*w^4-76*w^3+176*w^2)*x^6+(263-841/3*w+10*w^4-56*w^3+331/2*w^2)*x^5+(292-1385/6*w+105*w^2+4*w^4-30*w^3)*x^4+(821/4+w^4-733/6*w-10*w^3+91/2*w^2)*x^3+(173/2-2*w^3-122/3*w+11*w^2)*x^2+(77/4+3/2*w^2-37/6*w)*x+3/2-1/2*w)*3^(1/2)-1/4-3*w^4*x^12+(-12*w^4+10*w^3)*x^11+(-30*w^4+26*w^3-7*w^2)*x^10+(-48*w^4+54*w^3+14*w^2-14*w)*x^9+(63/2-57*w^4+64*w^3+36*w^2-98*w)*x^8+(140-48*w^4+68*w^3+80*w^2-565/3*w)*x^7+(255-769/3*w-30*w^4+48*w^3+127/2*w^2)*x^6+(292-1169/6*w-12*w^4+34*w^3+45*w^2)*x^5+(797/4-661/6*w+7/2*w^2-3*w^4+14*w^3)*x^4+(173/2-80/3*w+6*w^3-w^2)*x^3+(69/4-25/6*w-9/2*w^2)*x^2+(3/2+3/2*w)*x))/(x^2+x+1)^4/x^4;

(1/8)*2^(1/2)*(-2*x*(w^4*x^11+(4*w^4-6*w^3)*x^10+(10*w^4-22*w^3+17*w^2)*x^9+(16*w^4-50*w^3+62*w^2-30*w)*x^8+(71/2+19*w^4-72*w^3+132*w^2-114*w)*x^7+(140-(685/3)*w+16*w^4-76*w^3+176*w^2)*x^6+(263-(841/3)*w+10*w^4-56*w^3+(331/2)*w^2)*x^5+(292-(1385/6)*w+105*w^2+4*w^4-30*w^3)*x^4+(821/4+w^4-(733/6)*w-10*w^3+(91/2)*w^2)*x^3+(173/2-2*w^3-(122/3)*w+11*w^2)*x^2+(77/4+(3/2)*w^2-(37/6)*w)*x+3/2-(1/2)*w)*(-1+I*3^(1/2))^(1/2)+2^(1/2)*(I*x*(w^4*x^11+(4*w^4-6*w^3)*x^10+(10*w^4-22*w^3+17*w^2)*x^9+(16*w^4-50*w^3+62*w^2-30*w)*x^8+(71/2+19*w^4-72*w^3+132*w^2-114*w)*x^7+(140-(685/3)*w+16*w^4-76*w^3+176*w^2)*x^6+(263-(841/3)*w+10*w^4-56*w^3+(331/2)*w^2)*x^5+(292-(1385/6)*w+105*w^2+4*w^4-30*w^3)*x^4+(821/4+w^4-(733/6)*w-10*w^3+(91/2)*w^2)*x^3+(173/2-2*w^3-(122/3)*w+11*w^2)*x^2+(77/4+(3/2)*w^2-(37/6)*w)*x+3/2-(1/2)*w)*3^(1/2)-1/4-3*w^4*x^12+(-12*w^4+10*w^3)*x^11+(-30*w^4+26*w^3-7*w^2)*x^10+(-48*w^4+54*w^3+14*w^2-14*w)*x^9+(63/2-57*w^4+64*w^3+36*w^2-98*w)*x^8+(140-48*w^4+68*w^3+80*w^2-(565/3)*w)*x^7+(255-(769/3)*w-30*w^4+48*w^3+(127/2)*w^2)*x^6+(292-(1169/6)*w-12*w^4+34*w^3+45*w^2)*x^5+(797/4-(661/6)*w+(7/2)*w^2-3*w^4+14*w^3)*x^4+(173/2-(80/3)*w+6*w^3-w^2)*x^3+(69/4-(25/6)*w-(9/2)*w^2)*x^2+(3/2+(3/2)*w)*x))/((x^2+x+1)^4*x^4)

solve(eq = 0,w)

(1/2)*(2*x^2+1)/((x^2+x+1)*x)

PDEtools:-Solve(eq=0,w)

w = (1/2)*(2*x^2+1)/((x^2+x+1)*x)

restart;

interface(warnlevel=4);
kernelopts('assertlevel'=2):

3

eq:=1/8*2^(1/2)*(-2*x*(w^4*x^11+(4*w^4-6*w^3)*x^10+(10*w^4-22*w^3+17*w^2)*x^9+(16*w^4-50*w^3+62*w^2-30*w)*x^8+(71/2+19*w^4-72*w^3+132*w^2-114*w)*x^7+(140-685/3*w+16*w^4-76*w^3+176*w^2)*x^6+(263-841/3*w+10*w^4-56*w^3+331/2*w^2)*x^5+(292-1385/6*w+105*w^2+4*w^4-30*w^3)*x^4+(821/4+w^4-733/6*w-10*w^3+91/2*w^2)*x^3+(173/2-2*w^3-122/3*w+11*w^2)*x^2+(77/4+3/2*w^2-37/6*w)*x+3/2-1/2*w)*(-1+I*3^(1/2))^(1/2)+2^(1/2)*(I*x*(w^4*x^11+(4*w^4-6*w^3)*x^10+(10*w^4-22*w^3+17*w^2)*x^9+(16*w^4-50*w^3+62*w^2-30*w)*x^8+(71/2+19*w^4-72*w^3+132*w^2-114*w)*x^7+(140-685/3*w+16*w^4-76*w^3+176*w^2)*x^6+(263-841/3*w+10*w^4-56*w^3+331/2*w^2)*x^5+(292-1385/6*w+105*w^2+4*w^4-30*w^3)*x^4+(821/4+w^4-733/6*w-10*w^3+91/2*w^2)*x^3+(173/2-2*w^3-122/3*w+11*w^2)*x^2+(77/4+3/2*w^2-37/6*w)*x+3/2-1/2*w)*3^(1/2)-1/4-3*w^4*x^12+(-12*w^4+10*w^3)*x^11+(-30*w^4+26*w^3-7*w^2)*x^10+(-48*w^4+54*w^3+14*w^2-14*w)*x^9+(63/2-57*w^4+64*w^3+36*w^2-98*w)*x^8+(140-48*w^4+68*w^3+80*w^2-565/3*w)*x^7+(255-769/3*w-30*w^4+48*w^3+127/2*w^2)*x^6+(292-1169/6*w-12*w^4+34*w^3+45*w^2)*x^5+(797/4-661/6*w+7/2*w^2-3*w^4+14*w^3)*x^4+(173/2-80/3*w+6*w^3-w^2)*x^3+(69/4-25/6*w-9/2*w^2)*x^2+(3/2+3/2*w)*x))/(x^2+x+1)^4/x^4:

solve(eq = 0,w)

Error, (in Internal:-FactorEasy) assertion failed, Internal:-FactorEasy expects its return value to be of type list(_POWER(POLYNOMIAL, posint)), but computed [_POWER(POLYNOMIAL([0, [z, x, r1, r2, r3], [[[[1], [0, -1]], 0, [[1]]], [[1], 0, [1]], [-3, 0, 1]]], [[[0, [[1]]]]]), 1), _POWER(POLYNOMIAL([0, [z, x, r1, r2, r3], [[[[1], [0, -1]], 0, [[1]]], [[1], 0, [1]], [-3, 0, 1]]], [[[[0, [1]]]]]), FAIL), _POWER(POLYNOMIAL([0, [z, x, r1, r2, r3], [[[[1], [0, -1]], 0, [[1]]], [[1], 0, [1]], [-3, 0, 1]]], [[[[[0, 1]]]]]), FAIL), _POWER(POLYNOMIAL([0, [z, x, r1, r2, r3], [[[[1], [0, -1]], 0, [[1]]], [[1], 0, [1]], [-3, 0, 1]]], [[[[[-1], [0, 1]]], 0, [[[-8], [0, 8]]], 0, [[[-24], [0, 24]]], 0, [[[-32], [0, 32]]], 0, [[[-16], [0, 16 ... , [[[-160], [0, 160]]], [[[-64], [0, 64]]], [[[-16], [0, 16]]]]]), 1)]

PDEtools:-Solve(eq=0,w)

Error, (in Internal:-FactorEasy) assertion failed, Internal:-FactorEasy expects its return value to be of type list(_POWER(POLYNOMIAL, posint)), but computed [_POWER(POLYNOMIAL([0, [z, x, r1, r2, r3], [[[[1], [0, -1]], 0, [[1]]], [[1], 0, [1]], [-3, 0, 1]]], [[[0, [[1]]]]]), 1), _POWER(POLYNOMIAL([0, [z, x, r1, r2, r3], [[[[1], [0, -1]], 0, [[1]]], [[1], 0, [1]], [-3, 0, 1]]], [[[[0, [1]]]]]), FAIL), _POWER(POLYNOMIAL([0, [z, x, r1, r2, r3], [[[[1], [0, -1]], 0, [[1]]], [[1], 0, [1]], [-3, 0, 1]]], [[[[[0, 1]]]]]), FAIL), _POWER(POLYNOMIAL([0, [z, x, r1, r2, r3], [[[[1], [0, -1]], 0, [[1]]], [[1], 0, [1]], [-3, 0, 1]]], [[[[[-1], [0, 1]]], 0, [[[-8], [0, 8]]], 0, [[[-24], [0, 24]]], 0, [[[-32], [0, 32]]], 0, [[[-16], [0, 16 ... , [[[-160], [0, 160]]], [[[-64], [0, 64]]], [[[-16], [0, 16]]]]]), 1)]

 

Download assertion_failed.mw

restart:

sys:={-diff(v(x,t),t)+0.5*p*diff(u(x,t),x,x)+q*u(x,t)*(u(x,t)^2+v(x,t)^2)=0,diff(u(x,t),t)+0.5*p*diff(v(x,t),x,x)+q*v(x,t)*(u(x,t)^2+v(x,t)^2)=0};
        /                      /  2         \
        |/ d         \         | d          |
sys := < |--- u(x, t)| + 0.5 p |---- v(x, t)|
        |\ dt        /         |   2        |
        \                      \ dx         /

               /       2          2\       / d         \
   + q v(x, t) \u(x, t)  + v(x, t) / = 0, -|--- v(x, t)|
                                           \ dt        /

           /  2         \                                      \ 
           | d          |             /       2          2\    | 
   + 0.5 p |---- u(x, t)| + q u(x, t) \u(x, t)  + v(x, t) / = 0 >
           |   2        |                                      | 
           \ dx         /                                      / 
eq1 := diff(u(x,t),t) = u__t(x,t):
eq2 := diff(v(x,t),t) = v__t(x,t):

sys_tmp := subs(eq1, eq2, sys):

sys_new := sys_tmp union {eq1, eq2}:

Boundary conditions:
bc :=
    u(0,t) = 2,
    v(0,t) = 0;
                 bc := u(0, t) = 2, v(0, t) = 0
Initial conditions:
ic :=
    u(x,0) = tanh(2*Pi),
    v(x,0) = tanh(2*Pi),
    u__t(x,0) = 0,
    v__t(x,0) = 0;
 
ic := u(x, 0) = tanh(2 Pi), v(x, 0) = tanh(2 Pi), u__t(x, 0) = 0, 

  v__t(x, 0) = 0
Solve the system:
pdsol := pdsolve(subs(p=1, q=0.5, sys_new), {ic, bc}, numeric);

I want to change my variables. The coding below shows the influence of the Nb and Nt on the local Nusselt number (-theta') against the convection Bi.

But now I want it to be  the influence of the Bi on the local Nusselt number (-theta') against the Nb and Nt.

Can anyone help me?

NULL

with(plots)

DE1 := diff(f(eta), `$`(eta, 3))+f(eta)*(diff(f(eta), `$`(eta, 2)))-(diff(f(eta), eta))^2 = 0

diff(diff(diff(f(eta), eta), eta), eta)+f(eta)*(diff(diff(f(eta), eta), eta))-(diff(f(eta), eta))^2 = 0

(1)

DE2 := diff(theta(eta), `$`(eta, 2))+Pr*f(eta)*(diff(theta(eta), eta))+Pr*Nb*(diff(theta(eta), eta))*(diff(phi(eta), eta))+Pr*Nt*(diff(theta(eta), eta))^2 = 0

diff(diff(theta(eta), eta), eta)+Pr*f(eta)*(diff(theta(eta), eta))+Pr*Nb*(diff(theta(eta), eta))*(diff(phi(eta), eta))+Pr*Nt*(diff(theta(eta), eta))^2 = 0

(2)

DE3 := diff(phi(eta), `$`(eta, 2))+Le*(diff(phi(eta), eta))+Nt*(diff(theta(eta), `$`(eta, 2)))/Nb

diff(diff(phi(eta), eta), eta)+Le*(diff(phi(eta), eta))+Nt*(diff(diff(theta(eta), eta), eta))/Nb

(3)

BC1 := f(0) = 0, (D(f))(0) = 1, (D(f))(10) = 0

f(0) = 0, (D(f))(0) = 1, (D(f))(10) = 0

(4)

BC2 := theta(10) = 0, (D(theta))(0) = -Bi*(1-theta(0))

theta(10) = 0, (D(theta))(0) = -Bi*(1-theta(0))

(5)

BC3 := phi(0) = 1, phi(10) = 0

phi(0) = 1, phi(10) = 0

(6)

``

  getRes:= proc(a, x);
                local sol;
                if   type(x, numeric)
                then sol:= dsolve
                           ( eval
                             ( {BC1, BC2, BC3, DE1, DE2, DE3},
                               [Nt = a, Nb = a, Pr = 2, Le = 5, Bi = x]
                             ),
                             numeric,
                             output = listprocedure,
                             abserr = 0.0001,
                             maxmesh = 1024,
                             initmesh = 512
                           );
                else return 'procname(a,x)';
                fi;
                return -eval(diff(theta(eta), eta), sol)(0)
           end proc:

  cl:=[red, green, blue, yellow, black]:
  L:=[0.2, 0.4, 0.6,0.8,1.0]:
  display
  ( [ seq
      ( plot
        ( getRes( L[k], x),
          x=0..5,color = colorList[k], legend = ["Nb=Nt" = L[k]], legendstyle = [location = top],
          color=cl[k]
        ),
        k=1..5
      )
    ]
  )
               

 

 

 

Download odeProb2.mw

1 2 3 4 5 6 7 Last Page 1 of 1878