tomleslie

13771 Reputation

19 Badges

14 years, 13 days

MaplePrimes Activity


These are replies submitted by tomleslie

@AHSAN 

you will have with me is via this website

@AHSAN 

on my most recent post - or are you replying to yourself for some reason?

@AHSAN 

How quickly you would have received the answer you want, if only you had been able to ask the question clearly and unambiguously - because this is much simpler. See the attached

restart

n := 0

"u(y):=-(3 We (y-1-(x^2)/2) (y+1+(x^2)/2) (-(9 ((1+(x^2)/2)^2+y^2) (n-1) (k (1+(x^2)/2)+2 Q-1-(x^2)/2)^3)/(4 (1+(x^2)/2)^7)-(3 y (n-1) (k+1) (k (1+(x^2)/2)+2 Q-1-(x^2)/2)^2)/((1+(x^2)/2)^5)-(3 (k+1)^2 (n-1) (k (1+(x^2)/2)+2 Q-1-(x^2)/2))/(2 (1+(x^2)/2)^3)-(-120 lambda (1+(x^2)/2)^4-63 (k^2-4/7 k+1) (k-1) (n-1) (1+(x^2)/2)^3-288 Q (n-1) (k^2-11/8 k+1) (1+(x^2)/2)^2-486 Q^2 (n-1) (k-1) (1+(x^2)/2)-324 Q^3 (n-1))/(15 (1+(x^2)/2)^5)))/(16 (1+(x^2)/2)^2)+((3 k (1+(x^2)/2))/2+3 Q-3/2-(3 x^2)/4+(-(3 y^2 (k (1+(x^2)/2)+2 Q-1-(x^2)/2))/(2 (1+(x^2)/2)^3)-k+1) (1+(x^2)/2)-(k+1) y)/(2 (1+(x^2)/2)) :"

ode := diff(theta(y), y, y)+Br*((diff(u(y), y))^2+(1/2)*We*(n-1)*(diff(u(y), y))^4) = 0

bc := theta(-1-(1/2)*x^2) = 0, theta(1+(1/2)*x^2) = 1

Q__vals := [.5617, .4392, .2564, .1645, 0.659e-1]

`λ__vals` := [0.96e-2, 0.10e-2, 0.93e-2, 0.75e-2, 0.31e-2]; k__vals := [.1, .3, .5, .7, .9]

colors := [black, red, blue, green, orange, cyan, purple, yellow, navy]; for j to numelems(k__vals) do sol := dsolve(eval([ode, bc], [k = k__vals[j], We = .1, Br = .3, x = 0, Q = Q__vals[j], lambda = `λ__vals`[j]]), numeric); plt[j] := plots:-odeplot(sol, [y, theta(y)], y = -1 .. 1, axes = boxed, color = colors[j], legend = typeset("k=", k__vals[j])) end do; plots:-display([seq(plt[j], j = 1 .. 5)], title = typeset("We=%1, Br =%2, x =%3", .1, .3, 0))

 

``

Download again.mw

@AHSAN 

I have adjusted the coloring.

The orientation can be set using the orientation=[angle1, angle2, angle2] as an option to the display() command. One possibility is shown in the attached. If I am honest, working out which is the "best" list of angles, so the pragmatic way to do this is to left-click on the plot, hold the left mouse button down and just drag the cursor around until you get the view that you want. Then you can then read the orientation in the plot toolbar which is at the top of the screen. Insert these values in the orientation option, for teh display command. Any time you re-execute the worksheet, it will appear in the the orientation you specified

restart

psi := -(1.600000000*10^(-12))*(-(1.100464934*10^15)*x^6*y-(5.933879812*10^14)*x^8*y-(7.369485980*10^14)*x^2*y+(3.524964826*10^14)*x^2*y^2+(6.165015570*10^14)*x^4*y^2+(6.152509407*10^14)*x^6*y^2+(1.335184945*10^13)*x^2*y^3-(1.208534863*10^15)*x^4*y+(1.495218750*10^12)*x^16*y+(5.500000000*10^12)*x^14*y^2+(3.665811328*10^12)*x^14*y+(8.779561331*10^13)*y^2+(1.021933440*10^11)*y^4+(1.638534652*10^13)*y^3-(1.932126733*10^14)*y+(1.132673864*10^10)*y^5-(4.665845924*10^13)*x^4*y^3-(3.504346560*10^11)*x^2*y^4-(9.402963557*10^13)*x^6*y^3+(7.410783600*10^10)*x^4*y^4+(3.880305000*10^11)*x^6*y^4+(3.839362913*10^14)*x^8*y^2-(2.158659379*10^13)*x^12*y-(1.807479066*10^14)*x^10*y+(3.437500000*10^11)*x^16*y^2+(3.843735156*10^13)*x^12*y^2-(5.625000000*10^11)*x^14*y^3+(1.406250000*10^11)*x^18*y-(7.445478720*10^13)*x^8*y^3-(3.061199062*10^13)*x^10*y^3+(1.535553910*10^14)*x^10*y^2-(1.230187500*10^11)*x^6*y^5+(1.666494000*10^11)*x^4*y^5-(7.525146240*10^10)*x^2*y^5+(1.252968750*10^11)*x^8*y^4-(6.493625000*10^12)*x^12*y^3)/(x^2+2.)^9-(.2460375000*((-6.737387592+2.235431590*y)*x^16+(-3.657978966*y^2+35.76690545*y-107.8482015)*x^14+(-42.22852207*y^2+249.9609307*y-754.9196770)*x^12+(-199.0720316*y^2+998.5820269*y-3019.038361)*x^10+(.8148148150*y^3-484.1849698*y^2+2496.766005*y-7546.088582)*x^8+(-.8000000000*y^4+2.523390946*y^3-611.4816515*y^2+4001.022222*y-12072.72831)*x^6+(1.083733335*y^4+.4819287219*y^3-303.4233996*y^2+4009.155072*y-12073.18307)*x^4+(-.4893657917*y^4-2.278902388*y^3+86.82806144*y^2+2292.310612*y-6899.426943)*x^2+0.7365863430e-1*y^4+.6645708495*y^3+106.5551164*y^2+570.9413454*y-1724.697565))*y/(x^2+2.)^9

-0.1600000000e-11*(-0.1932126733e15*y+0.8779561331e14*y^2+0.1021933440e12*y^4+0.1638534652e14*y^3+0.1132673864e11*y^5+0.1666494000e12*x^4*y^5-0.7525146240e11*x^2*y^5-0.5625000000e12*x^14*y^3-0.6493625000e13*x^12*y^3-0.3061199062e14*x^10*y^3+0.1252968750e12*x^8*y^4-0.1230187500e12*x^6*y^5+0.1406250000e12*x^18*y+0.3437500000e12*x^16*y^2+0.3843735156e14*x^12*y^2-0.2158659379e14*x^12*y+0.1535553910e15*x^10*y^2-0.1807479066e15*x^10*y-0.7445478720e14*x^8*y^3+0.3839362913e15*x^8*y^2+0.3880305000e12*x^6*y^4-0.9402963557e14*x^6*y^3+0.7410783600e11*x^4*y^4-0.4665845924e14*x^4*y^3-0.3504346560e12*x^2*y^4+0.1335184945e14*x^2*y^3+0.1495218750e13*x^16*y+0.5500000000e13*x^14*y^2+0.3665811328e13*x^14*y+0.6152509407e15*x^6*y^2+0.6165015570e15*x^4*y^2+0.3524964826e15*x^2*y^2-0.5933879812e15*x^8*y-0.1100464934e16*x^6*y-0.1208534863e16*x^4*y-0.7369485980e15*x^2*y)/(x^2+2.)^9-.2460375000*((-6.737387592+2.235431590*y)*x^16+(-3.657978966*y^2+35.76690545*y-107.8482015)*x^14+(-42.22852207*y^2+249.9609307*y-754.9196770)*x^12+(-199.0720316*y^2+998.5820269*y-3019.038361)*x^10+(.8148148150*y^3-484.1849698*y^2+2496.766005*y-7546.088582)*x^8+(-.8000000000*y^4+2.523390946*y^3-611.4816515*y^2+4001.022222*y-12072.72831)*x^6+(1.083733335*y^4+.4819287219*y^3-303.4233996*y^2+4009.155072*y-12073.18307)*x^4+(-.4893657917*y^4-2.278902388*y^3+86.82806144*y^2+2292.310612*y-6899.426943)*x^2+0.7365863430e-1*y^4+.6645708495*y^3+106.5551164*y^2+570.9413454*y-1724.697565)*y/(x^2+2.)^9

(1)

with(plots)

p1 := contourplot(psi, x = -3 .. 3, y = -.5 .. .5, coloring = ["Blue", "Grey"], axes = boxed, filledregions = true)

 

f := plottools:-transform(proc (x, y) options operator, arrow; [x, y, -1] end proc)

display([contourplot3d(psi, x = -3 .. 3, y = -.5 .. .5, coloring = ["Blue", "Grey"], filledregions = true), f(p1)], orientation = [60, 60, 30])

 

NULL

NULL

Download contPlot3.mw

 

@AHSAN 

The attached still produces plots 25 at a time, but now there are only five curves on each plot. It akso shows how to display individual plots, depending on which values of 'Q' and 'lambda' you are interested in

restart

n := 0

"u(y):=-(3 We (y-1-(x^2)/2) (y+1+(x^2)/2) (-(9 ((1+(x^2)/2)^2+y^2) (n-1) (k (1+(x^2)/2)+2 Q-1-(x^2)/2)^3)/(4 (1+(x^2)/2)^7)-(3 y (n-1) (k+1) (k (1+(x^2)/2)+2 Q-1-(x^2)/2)^2)/((1+(x^2)/2)^5)-(3 (k+1)^2 (n-1) (k (1+(x^2)/2)+2 Q-1-(x^2)/2))/(2 (1+(x^2)/2)^3)-(-120 lambda (1+(x^2)/2)^4-63 (k^2-4/7 k+1) (k-1) (n-1) (1+(x^2)/2)^3-288 Q (n-1) (k^2-11/8 k+1) (1+(x^2)/2)^2-486 Q^2 (n-1) (k-1) (1+(x^2)/2)-324 Q^3 (n-1))/(15 (1+(x^2)/2)^5)))/(16 (1+(x^2)/2)^2)+((3 k (1+(x^2)/2))/2+3 Q-3/2-(3 x^2)/4+(-(3 y^2 (k (1+(x^2)/2)+2 Q-1-(x^2)/2))/(2 (1+(x^2)/2)^3)-k+1) (1+(x^2)/2)-(k+1) y)/(2 (1+(x^2)/2)) :"

ode := diff(theta(y), y, y)+Br*((diff(u(y), y))^2+(1/2)*We*(n-1)*(diff(u(y), y))^4) = 0

bc := theta(-1-(1/2)*x^2) = 0, theta(1+(1/2)*x^2) = 1

Q__vals := [.5617, .4392, .2564, .1645, 0.659e-1]

`λ__vals` := [0.96e-2, 0.10e-2, 0.93e-2, 0.75e-2, 0.31e-2]; arrPlt := Array(1 .. numelems(Q__vals), 1 .. numelems(`λ__vals`))

colors := [black, red, blue, green, orange, cyan, purple, yellow, navy]; for p to numelems(Q__vals) do for q to numelems(`λ__vals`) do for j by 2 to 9 do sol := dsolve(eval([ode, bc], [k = .1*j, We = .1, Br = .3, x = 0, Q = Q__vals[p], lambda = `λ__vals`[q]]), numeric); plt[j] := plots:-odeplot(sol, [y, theta(y)], y = -1 .. 1, axes = boxed, color = colors[j], legend = typeset("k=", .1*j)) end do; arrPlt[p, q] := plots:-display([seq(plt[j], j = 1 .. 9, 2)], title = typeset("We=%1, Br =%2, x =%3 Q=%4 λ=%5", .1, .3, 0, Q__vals[p], `λ__vals`[q])) end do end do

plots:-display(arrPlt[4, 1]); plots:-display(arrPlt[2, 3])

 

 

 

Download pltArr2.mw

@AHSAN 

You state

Can you please explain about Arrplt?

Arrplt is just a 2-dimensionl array which is used to store plots

You state

please convert this command into noraml display beacuse when i copy to paste into file its looks stratched.

Arrplt is not a command and I have no idea what "its looks stratched." means

You state

Sir, I observed that  ploted value of Q and Lmabda are againt for k=0.1,0.3,0.5,0.7,0.9 and We=0.1,0.3,0.5,0.7,0.9 but in graph you draw for k=0.1,0.2,,,0.9.

No! As I stated in my original response.

The code in the attached, uses

Q=[0.5617, 0.4392, 0.2564, 0.1645, 0.0659]
lambda=[0.0096, 0.0010, 0.0093, 0.0075, 0.0031]

to produce an array of 25 plots, organised as 5x5. Each individual plot contains nine curves, given by

k = 0.1*j, We = 0.1, Br = 0.3, x = 0  with j=1..9

@Zeineb 

The earliest version of Maple which runs Rouben's code without producing errors is Maple2018. All earlier versions ie (Maple2017, Maple2016, Maple2015, and Maple18) generate errors. Since you are using Maple 17 (ie an even earlier version, released in 2013), I think you can assume that it will generate errors as well.

That's a problem with using seriously out-of-date software

@AHSAN 

(amongst many).

  1. The left-hand-side of the first equation in your pdf does not correspond to the quivalent term in yur worksheet. Which i correct?
  2. In both the pdf and your worksheet, the LHS of the first equation "appears" to be a function of 'y' only, whereas the RHS is a function of 'x' only. How can that possibly be true if 'x' and 'y' are independent variables?
  3. I suggest you re-read bullet point (4) in my earlier response, because with n=2, it is impossible to generate a power series in 'we'

I don't see anything wrong with your approach in

sol:=isolate(-2*y-cos(y)+x+sin(x)-_C1 = 0,y);
if lhs(sol)=y then
   print("it solved it. Soltion is ",sol);
else
   print("Failed to solve it. Try solve command now ");
fi;

but I'd probably want to handle the case where the eqquation happened to be of the form y=f(x,y) - where 'y' is already on the LHS. For any equation to be tested, I'd probably convert it to the form lhs(eqn)-rhs(eqn)=0, before trying isolate()

@JAMET 

Is this a question? What problem do you have?

Quote from my previous answer

In future please use the big, green up-arrow in the Mapleprimes toolbar to upload a actual worksheet. When you cut-and-paste as  in your original question, responders have too go through it line-by-line to work out what is Maple inpu, what is Maple output, what is comment, what is plaintext etc etc.  This is very tedious!

@acer 
interesting!

  1. I didn't check the Physics Updates version for the first (Maple2022.1) worksheet which I posted, and I'm a bit "careless" about routinely installing the former. I *could* have been using Physics Updates from anything up to a couple of months ago. So, for me,  the backticks approach works in Maple 2022.1 with an *unknown* version of the Physics toolbox
  2. Maple 2022.2, which I installed today (shortly after midnight, 02/11/2022) produced the second worksheet I posted. Maaple2022.2 comes with a "new" set of Physics Updates, which appears to be version 1337 created 2022, October 19, 14:41 hours Pacific Time. Checking the "Date Modified" field in the appropriate toolbox directory confirms that a couple of these toolbox files were modified today (02/11/2022). So, for me,  the backticks approach works in Maple 2022.2 with Physics Updates version 1337
  3. Just for completeness, I have now installed the "latest and greatest" version of Physics Updates (1340) and the backticks approach is still working - see attached

  restart;
  interface(version);
  Physics:-Version();

`Standard Worksheet Interface, Maple 2022.2, Windows 7, October 23 2022 Build ID 1657361`

 

`The "Physics Updates" version in the MapleCloud is 1340 and is the same as the version installed in this computer, created 2022, October 31, 16:51 hours Pacific Time.`

(1)

#
# With 1-D input, use backticks
# Without them this fails
#
  X__a,b:=1;

Error, mismatched multiple assignment of 2 variables on the left side and 1 value on the right side

 

#
# But this does not
#
  `Y__a,b`:=1;
  `Z__a,b,c`:=1;

1

 

1

(2)

`X__a,b` := 1; `Y__a,b` := 1

1

 

1

(3)

NULL

Download backT3.mw

@acer 

I missed the differrence between 2022.1 and 2022.2 - but uing bacticks on input stoill seems to provide "correct" output.

See attached.

  restart;
  interface(version);

`Standard Worksheet Interface, Maple 2022.2, Windows 7, October 23 2022 Build ID 1657361`

(1)

#
# With 1-D input, use backticks
# Without them this fails
#
  X__a,b:=1;

Error, mismatched multiple assignment of 2 variables on the left side and 1 value on the right side

 

#
# But this does not
#
  `Y__a,b`:=1;
  `Z__a,b,c`:=1;

1

 

1

(2)

`X__a,b` := 1; `Y__a,b` := 1

1

 

1

(3)

NULL

Download backT2.mw

@romanrieme 

which does produce the animation. Two points

  1. I haven't even looked at the 'maxfun' problem, but this means that The the animation range is restricted to t<~0.7
  2. technique to produce the animation is a bit crude, but is governed by the way earlier code is written (and I haven't got the time for a major rewrite!)

restart

q1 := 5

b1 := 18

l1 := 16

r1 := .35

ms := 626

g1 := -9.807

u0 := -19.4

v0 := -17.7

 

x0 := 0

 

z0 := 15

 

xa := 0

 

za := z0+1

 

mu1 := .35

omega0 := -24

mn := b1*g1*q1*za+g1*mp

fz := mn/za

Js := (2/5)*ms*r1^2

"Lp(t):=sqrt((xn(t)-xa)^(2)+(zn(t)-za)^(2)): #"

"Lr(t):=piecewise((za-Lp(t))>0,za-Lp(t),0) :#"

"A1(t):=piecewise((xa-xn(t))>0,xa-xn(t),0) :#"

"sinalpha1(t):=piecewise((xa-xn(t))>0,(xa-xn(t))/(Lp(t)),0) :#"

"cosalpha1(t):=piecewise((xa-xn(t))>0,((za-zn(t)))/(Lp(t)),1) :#"

"alpha2(t):=arccos(cosalpha1(t))+phi(t) :#"

``

"S1(t):=piecewise(zn(t)>=0 and xn(t)<=xa,fz*Lr(t),0):#"

"S1z(t):=S1(t)*cos(phi(t)):#"

"S1x(t):=S1(t)*sin(-phi(t)):#"

"S2(t):=piecewise(xn(t)<=xa and zn(t)>0, (-1) S1(t)*exp(mu1*alpha2(t)),0):#"

"S2x(t):=S2(t)*sinalpha1(t):#"

"S2z(t):=S2(t)*cosalpha1(t):#"

"Sx(t):=S2x(t)+S1x(t):#"

"Sz(t):=S1z(t)+S2z(t)+ms*g1:#"

shape_factor := 0

"Fr(t):=S2(t)+S1(t)+shape_factor:"

NULL

"m(t):=Lr(t)q1*b1:#"

mp := 300

"Jp(t):=1/(3)*m(t)*Lr(t)^(2)+mp*Lr(t)^(2) :#"

"M(t):=g1*m(t)*(Lr(t))/(2)*sin(phi(t))+g1*mp*sin(phi(t)):#"

k := 800*ms

NULL

"Es(t):=(ms*g1*(z0-zn(t)))+1/(2)*ms*(vn(t)^(2)+un(t)^(2))+1/(2)*Js*omega(t)^(2):"

"En(t):=(1)/(2)*Jp(t)*((&DifferentialD;)/(&DifferentialD; t)phi(t))^(2)+-(m(t)+mp)*g1*(Lr(t)-cos(phi(t))*Lr(t)):"

NULL

eqx1 := ms*(diff(u(t), t)) = 0.; eqz1 := ms*(diff(v(t), t)) = ms*g1

eqx0 := diff(x(t), t) = u(t); eqz0 := diff(z(t), t) = v(t)

``

eqx1n := ms*(diff(un(t), t)) = Sx(t); eqz1n := ms*(diff(vn(t), t)) = ms*g1+S1z(t)+S2z(t)

eqx0n := diff(xn(t), t) = un(t); eqz0n := diff(zn(t), t) = vn(t)

eqr1 := Js*(diff(omega(t), t)) = (S2(t)+S1(t)+shape_factor)*r1

u0n := u0*ms/(ms+(1/3)*m(t)+mp)

u0n1 := u0*ms/(ms+mn)

phid00 := u0n1/(l1-z0)

phid0 := .7*eval[recurse](-u0n(t)/(l1-z0), [t = 0, xn(0) = x0, zn(0) = z0])

simplify(phid0)

"znb(t):=(Lr(t)-(sin((Pi)/(2)- phi(t))*Lr(t))):"

"xnb(t):=Lr(t)*cos(Pi/(2)-phi(t)):"

eqp := (diff(phi(t), t, t))*Jp(t) = M(t)-k*(diff(phi(t), t))

with(plots)

ini := x(0) = x0, u(0) = u0, z(0) = z0, v(0) = v0, xn(0) = x0, un(0) = u0, zn(0) = z0, vn(0) = v0, omega(0) = omega0, phi(0) = 0, (D(phi))(0) = phid0

sol1 := dsolve({eqp, eqr1, eqx0, eqx0n, eqx1, eqx1n, eqz0, eqz0n, eqz1, eqz1n, ini}, {omega(t), phi(t), u(t), un(t), v(t), vn(t), x(t), xn(t), z(t), zn(t)}, type = numeric, output = operator)

NULL

odeplot(sol1, [[xn(t), zn(t)]], 0 .. 2.5, titlefont = ["HELVETICA", bolditalic, 10], labels = ["x (m)", "z (m)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], linestyle = [dashdot, solid], size = [300, 400])

Warning, cannot evaluate the solution further right of .70357836, maxfun limit exceeded (see ?dsolve,maxfun for details)

 

 

odeplot(sol1, [xnb(t), znb(t)], 0 .. 1, titlefont = ["HELVETICA", bolditalic, 10], size = [400, 400])

Warning, cannot evaluate the solution further right of .70357836, maxfun limit exceeded (see ?dsolve,maxfun for details)

 

 

with(plottools)

`~`[assign](sol1)

pt1 := [xa, za]

"pt2(t):=[xn(t),zn(t)]:"

"pt3(t):=[xnb(t), znb(t)]:"

animate(display, [line(pt1, pt2(t), color = blue, thickness = 4), line(pt2(t), pt3(t), color = red, thickness = 4), 'plot([pt2(t)], style = point, symbol = solidcircle, symbolsize = 16, color = blue)', 'plot([pt3(t)], style = point, symbol = solidcircle, symbolsize = 16, color = red)'], t = 0 .. .7, background = pointplot([pt1], symbol = solidcircle, symbolsize = 24, color = green))

 

``

Download anim2.mw

@acer 

test2() calls test1() with arguments 2, 3=3. so test1 computes 2+1 and returns 3, so test2 returns 2*3, ie 6

@acer 

the naming 'clashes' and their solution - but I am really curious how the OP's original worksheet manged to return the value 4, for the invocations

test2(2, b=3);

and

test3(2, b=3);

because I would have placed a large bet on a retrun value of 6, in both cases.

Indeed, I have just run the OP's original code and I do get a return value of 6 in both cases. See the attached.

So how did (s)he manage to get 4 ?????

restart:
interface(version);

`Standard Worksheet Interface, Maple 2022.1, Windows 7, May 26 2022 Build ID 1619613`

(1)

test1 := proc( a :: posint, { b :: posint := 1 } ) :: posint:
               return( a + b ):
         end proc:
test2 := proc( a :: posint, { b :: posint := 1 } ) :: posint:
               return( a * test1( a, b = b ) ):
         end proc:

test2(2, b=3);

6

(2)

test1(2, b=3);

5

(3)

2 * test1(2, b=3);

10

(4)

test3 := proc( a :: posint, { b :: posint := 1 } ) :: posint:
               return( a * test1( a, `b` = b ) ):
         end proc:

test3(2, b=3);

6

(5)

test4 := proc( a :: posint, { c :: posint := 1 } ) :: posint:
        return( a * test1( a, b = c ) ):
end proc:

test4(2, c=3);

10

(6)

 

 

Download keyParam.mw

First 9 10 11 12 13 14 15 Last Page 11 of 207