Advanced Search

tomleslie

13428 Reputation

19 Badges

13 years, 41 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

You are using...

tomleslie 13428
March 22 2023
0 1

the procedures f1() and f2() in generating equations in e1.. e6, however you are calling these with "square" brackets, eg f1[n]. These should be "round" brackets. In Maple "square" brackets are only used to access entries in an indexable data container

I have fixed this problem in the attached, where I have also restricted the extent of the main "for" loop, just to avoid generating too much output for this page,. Also, just for display purposes, I slightly modified the format strings in the main printf() statement, just for clarity.

For some reason this worksheet won't display inline here, so will you have to download to check

fsol.mw

Maybe...

tomleslie 13428
March 22 2023
1 0

using PDEtools:-dchange(), as shown in the attached?

> 

restart
> 

"T(x,eta):=Te+Te*(1-g(x,eta))*phi(x); PDE:=diff(T(x,eta),eta$2)+ f(x,eta)*diff(T(x,eta),eta)=x*(diff(f(x,eta),eta)*diff(T(x,eta),x)-diff(T(x,eta),eta)*diff(f(x,eta),x)) "

proc (x, eta) options operator, arrow, function_assign; Te+Te*(1-g(x, eta))*phi(x) end proc

  

-Te*(diff(diff(g(x, eta), eta), eta))*phi(x)-f(x, eta)*Te*(diff(g(x, eta), eta))*phi(x) = x*((diff(f(x, eta), eta))*(-Te*(diff(g(x, eta), x))*phi(x)+Te*(1-g(x, eta))*(diff(phi(x), x)))+Te*(diff(g(x, eta), eta))*phi(x)*(diff(f(x, eta), x)))

 (1)
> 

> 

NULL
> 

NULL
> 

NULL
> 

PDEtools:-dchange(eta = (ue/(c*x))^(1/2)*y, PDE, [y], params = [ue, c])

-Te*(diff(diff(g(y, x), y), y))*c*x*phi(x)/ue-f(y, x)*Te*(diff(g(y, x), y))*phi(x)/(ue/(c*x))^(1/2) = x*((diff(f(y, x), y))*(-Te*(diff(g(y, x), x))*phi(x)+Te*(1-g(y, x))*(diff(phi(x), x)))/(ue/(c*x))^(1/2)+Te*(diff(g(y, x), y))*phi(x)*(diff(f(y, x), x))/(ue/(c*x))^(1/2))

 (2)

NULL

Download chngeVar.mw

Or...

tomleslie 13428
March 21 2023
0 0

you could use complexplot(), as in the attached

> 

restart
> 

with(plottools)
> 

with(plots)
> 

with(CurveFitting)
> 

Digits := 10
> 

with(GaussInt)
> 

w := GInearest(0+I)

I

 (1)
> 

NULL
> 

"f(t):=7.0*(e)^((-(t-13180)^(2))/(2000000))+4.7*(e)^((-(t-16000)^(2))/(3200000)):"
> 

p1 := plot(f(t), t = 0 .. 20000, color = green); plots[display]({p1})

 
> 

NULL
> 

D1 := 15
> 

epsilon := 200000
> 

L := 6500
> 

v := .7

.7

 (2)
> 

t := 10000

10000

 (3)
> 

i = sqrt(-1)

i = I

 (4)
> 

"k(n) := evalf((2 *Pi*n)/(L))"

proc (n) options operator, arrow, function_assign; evalf(2*Pi*n/L) end proc

 (5)
> 

f(n) = (int(f(t)*exp(-w*k(n)*x), x = 0 .. L))/L
> 

``
> 

"C(x, t) := (∑) exp(-v* t*k(n)- D1 *t*(k(n)^())^(2)- epsilon *t*(k(n))^(4)) *f(n)* exp(w*k(n)* x)"

proc (x, t) options operator, arrow, function_assign; sum(exp(-v*t*k(n)-D1*t*k(n)^2-varepsilon*t*k(n)^4)*f(n)*exp(w*k(n)*x), n = 1 .. 10) end proc

 (6)
> 

uu10000 := [seq(evalf(C(L-j, t)), j = 0 .. 6500, 100)]; complexplot(uu10000, x = min(`~`[Re](uu10000)) .. max(`~`[Re](uu10000)), scaling = constrained, style = point)

 
> 

``

Download cplxplt.mw

The OP's primary objective...

tomleslie 13428
March 19 2023
0 0

is to solve the ODE given as f'''+ff''-f'^2-Mf'=0, with f(0)=0, f'(0)=1, and f'(5)=0 with M=0.5

So let's just solve it, as shown in the attached.

> 

  restart:
  with(plots):
> 

  M:=0.5;
  ode:= diff(f(x), x$3)+f(x)*diff(f(x), x$2)-diff(f(x),x)^2-M*diff(f(x),x)=0;
  conds:= f(0)=0, D(f)(0)=1, D(f)(5 )-0;
  ans:=dsolve( [ode, conds], numeric);
  odeplot(ans, [x, f(x)], x=0..5);
  

.5

  

diff(diff(diff(f(x), x), x), x)+f(x)*(diff(diff(f(x), x), x))-(diff(f(x), x))^2-.5*(diff(f(x), x)) = 0

  

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

  

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(13, {(1) = .0, (2) = .3648083980208248, (3) = .739453123146836, (4) = 1.1294231031240023, (5) = 1.5404396951722952, (6) = 1.9674194416734547, (7) = 2.4067451277606597, (8) = 2.850771504517652, (9) = 3.2971678108063687, (10) = 3.7442923732839364, (11) = 4.191663594165653, (12) = 4.620140941040035, (13) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(13, 3, {(1, 1) = .0, (1, 2) = 1.0, (1, 3) = -1.2249064035636923, (2, 1) = .294194453432655, (2, 2) = .6396128765313516, (2, 3) = -.7836115514965457, (3, 1) = .48635488180500724, (3, 2) = .40415085049570865, (3, 3) = -.49534834686074514, (4, 1) = .6116367675260377, (4, 2) = .250541601628673, (4, 3) = -.3073631245893844, (5, 1) = .6925060924519302, (5, 2) = .15125045647457197, (5, 3) = -.18593456512252354, (6, 1) = .7427378377645504, (6, 2) = 0.893905555320538e-1, (6, 3) = -.110380533553483, (7, 1) = .7730294979143654, (7, 2) = 0.51843487043024994e-1, (7, 3) = -0.6463752743281678e-1, (8, 1) = .7906832855559992, (8, 2) = 0.29654107059807183e-1, (8, 3) = -0.3774106284132804e-1, (9, 1) = .8007523087394047, (9, 2) = 0.1661644456196504e-1, (9, 3) = -0.22098378262475696e-1, (10, 1) = .8063158313666053, (10, 2) = 0.8938405451087109e-2, (10, 3) = -0.1307531449846363e-1, (11, 1) = .8092030889605896, (11, 2) = 0.4352616339821983e-2, (11, 3) = -0.7907830531921609e-2, (12, 1) = .810440946160643, (12, 2) = 0.1626507257535621e-2, (12, 3) = -0.5077258781403503e-2, (13, 1) = .8107323262402987, (13, 2) = .0, (13, 3) = -0.3612766608982084e-2}, datatype = float[8], order = C_order); YP := Matrix(13, 3, {(1, 1) = 1.0, (1, 2) = -1.2249064035636923, (1, 3) = 1.5, (2, 1) = .6396128765313516, (2, 2) = -.7836115514965457, (2, 3) = .9594452421864268, (3, 1) = .40415085049570865, (3, 2) = -.49534834686074514, (3, 3) = .6063284218940224, (4, 1) = .250541601628673, (4, 2) = -.3073631245893844, (4, 3) = .37603648294155106, (5, 1) = .15125045647457197, (5, 2) = -.18593456512252354, (5, 3) = .2272627479658001, (6, 1) = 0.893905555320538e-1, (6, 2) = -.110380533553483, (6, 3) = .13466974800716747, (7, 1) = 0.51843487043024994e-1, (7, 2) = -0.6463752743281678e-1, (7, 3) = 0.7857620604810917e-1, (8, 1) = 0.29654107059807183e-1, (8, 2) = -0.3774106284132804e-1, (8, 3) = 0.4554764716317479e-1, (9, 1) = 0.1661644456196504e-1, (9, 2) = -0.22098378262475696e-1, (9, 3) = 0.26279655923937467e-1, (10, 1) = 0.8938405451087109e-2, (10, 2) = -0.1307531449846363e-1, (10, 3) = 0.1509193089776011e-1, (11, 1) = 0.4352616339821983e-2, (11, 2) = -0.7907830531921609e-2, (11, 3) = 0.8594294332320503e-2, (12, 1) = 0.1626507257535621e-2, (12, 2) = -0.5077258781403503e-2, (12, 3) = 0.4930717565329715e-2, (13, 1) = .0, (13, 2) = -0.3612766608982084e-2, (13, 3) = 0.29289866770633205e-2}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(13, {(1) = .0, (2) = .3648083980208248, (3) = .739453123146836, (4) = 1.1294231031240023, (5) = 1.5404396951722952, (6) = 1.9674194416734547, (7) = 2.4067451277606597, (8) = 2.850771504517652, (9) = 3.2971678108063687, (10) = 3.7442923732839364, (11) = 4.191663594165653, (12) = 4.620140941040035, (13) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(13, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = 0.1273975608644069e-9, (2, 1) = -0.9699958654485496e-7, (2, 2) = 0.11898009711795105e-6, (2, 3) = -0.14531470840848373e-6, (3, 1) = -0.10320306463506545e-6, (3, 2) = 0.12691811465663893e-6, (3, 3) = -0.154678061248875e-6, (4, 1) = -0.7871067918956964e-7, (4, 2) = 0.9744291713856059e-7, (4, 3) = -0.11818172410834043e-6, (5, 1) = -0.4865437588049475e-7, (5, 2) = 0.613589205139654e-7, (5, 3) = -0.7358044074399424e-7, (6, 1) = -0.21030326071357806e-7, (6, 2) = 0.28453569331208245e-7, (6, 3) = -0.3288511507213022e-7, (7, 1) = -0.12272437930054753e-8, (7, 2) = 0.5323548042522246e-8, (7, 3) = -0.4173055146179112e-8, (8, 1) = 0.9931557233651253e-8, (8, 2) = -0.7035041293820884e-8, (8, 3) = 0.11351363226824974e-7, (9, 1) = 0.14250962440952772e-7, (9, 2) = -0.10817614549677532e-7, (9, 3) = 0.1640579916342742e-7, (10, 1) = 0.14719058046583237e-7, (10, 2) = -0.965367179587185e-8, (10, 3) = 0.15469412664652983e-7, (11, 1) = 0.13753746410766736e-7, (11, 2) = -0.6455085789179595e-8, (11, 3) = 0.12146772829413411e-7, (12, 1) = 0.1268333027160663e-7, (12, 2) = -0.2912324487441674e-8, (12, 3) = 0.8508819311976027e-8, (13, 1) = 0.12136283725856026e-7, (13, 2) = .0, (13, 3) = 0.5688098437791177e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[13] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.54678061248875e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 13, [f(x), diff(f(x), x), diff(diff(f(x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[13] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[13] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(3, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(13, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(13, 3, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[f(x), diff(f(x), x), diff(diff(f(x), x), x)]'[i] = yout[i], i = 1 .. 3)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[13] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.54678061248875e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 13, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[13] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[13] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(13, 3, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(3, {(1) = 0., (2) = 0., (3) = 0.}); `dsolve/numeric/hermite`(13, 3, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 3)] end proc, (2) = Array(0..0, {}), (3) = [x, f(x), diff(f(x), x), diff(diff(f(x), x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[f(x), diff(f(x), x), diff(diff(f(x), x), x)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc

  

 
> 

 

Download solODE.mw

8-th order polynomial...

tomleslie 13428
March 17 2023
0 0

which if you differentiate it, will give a 7-th order polynomial. Setting the latter equal to zero will rsult in seven roots (ie seven stationary points), although as mmcdara has pointed out three of these occur at x=-1. \Whether you consider this a single stationary point (with a multiplicity of 3), or three coincident stationary points is a quibble about language

The other four roots are complex

Second vertical axis is not logarithmic...

tomleslie 13428
March 17 2023
1 1

Check the tickmarks!

The problem appears to be that you have two sets of horizontal gridlines - one for each of the vertical axes in the dualaxis plot. Since the left and right tickmarks do not align, these two sets of gridlines are interleaved, giving an "illusion" of a weird vertical scale.

In a dualaxis plot, I can't think of any way to make the left and right tickmarks "align", since they are likely to be on very different scales.

Even if one could come up with a method of having horizontal gridlines only for the left axis (and I don't know if this can be done!), it means that these gridlines would not coincide with the tickmarks on the right axis - which would make the plot look a bit weird

anything...

tomleslie 13428
March 16 2023
0 0

wrong with the fsolve() command using the 'complex' option?

You couldread this response for some sug...

tomleslie 13428
March 14 2023
0 0

https://www.mapleprimes.com/questions/235892-How-Do-I-Solve-A-Differential-Equation-In-Maple

You may...

tomleslie 13428
March 14 2023
2 0

be able to solve it numerically, provided

  1. you supply numeric values for the parameters q1, q2, q3, q4
  2. you supply three initial/boundary conditions for the ODE

In the attached, I just use some random values for thes quantities

ode := -(q2-q4)*(q2-q3)*(q1-q4)*(q1-q3)*y(w)^2+(q1+q2-q3-q4)*(2*q1*q2-q1*q3-q1*q4-q2*q3-q2*q4+2*q3*q4)*(diff(y(w), w))*y(w)+(-2*q1^2-2*q1*q2+3*q1*q3+3*q1*q4-2*q2^2+3*q2*q3+3*q2*q4-6*q3*q4)*(diff(y(w), w))^2+(q1-q2+q3-q4)*(q1-q2-q3+q4)*(diff(y(w), w, w))*y(w)+(3*q1+3*q2-3*q3-3*q4)*(diff(y(w), w, w))*(diff(y(w), w))-3*(diff(y(w), w, w))^2+(-q1-q2+q3+q4)*(diff(y(w), w, w, w))*y(w)+2*(diff(y(w), w, w, w))*(diff(y(w), w)) = 0

-(q2-q4)*(q2-q3)*(q1-q4)*(q1-q3)*y(w)^2+(q1+q2-q3-q4)*(2*q1*q2-q1*q3-q1*q4-q2*q3-q2*q4+2*q3*q4)*(diff(y(w), w))*y(w)+(-2*q1^2-2*q1*q2+3*q1*q3+3*q1*q4-2*q2^2+3*q2*q3+3*q2*q4-6*q3*q4)*(diff(y(w), w))^2+(q1-q2+q3-q4)*(q1-q2-q3+q4)*(diff(diff(y(w), w), w))*y(w)+(3*q1+3*q2-3*q3-3*q4)*(diff(diff(y(w), w), w))*(diff(y(w), w))-3*(diff(diff(y(w), w), w))^2+(-q1-q2+q3+q4)*(diff(diff(diff(y(w), w), w), w))*y(w)+2*(diff(diff(diff(y(w), w), w), w))*(diff(y(w), w)) = 0

 (1)

params := [q1 = 1, q2 = -2, q3 = -3, q4 = 4]; ics := y(0) = 0, (D(y))(0) = 1, (D[1, 1](y))(0) = 0; ans := dsolve(eval([ode, ics], params), numeric); with(plots); odeplot(ans, [w, y(w)], w = 0 .. 5)

  

NULL

NULL

Download numODE.mw

Well...

tomleslie 13428
March 09 2023
0 0

You could start by typing

?How do I solve an ordinary differential equation

at the Maple prompt. If yoou need more info, check the Maple help for the command dsolve()

Once you have assigned...

tomleslie 13428
March 09 2023
1 0

values for 'm' and 'l' (I just used 1 for both), then you could try DEplot3d() as shown in the attached. I just guessed a value for the range of the dependent variable p(t).

 

For some reason, this particular plot won't render on this site, by I can assure you that this code

VF.mw

  restart;
  with(DETools):
  sys := { diff(r(t),t)=p(t)/m,
           diff(p(t),t)=l^2/(m*r(t)^3)-n*k*r(t)^(n-1),
           diff(phi(t),t)=l/(m*r(t)^2)
         };
  sys1:=subs({n=1,k=1, m=1, l=1},sys);
  p1:=DEplot3d( sys1,
                [r(t), p(t), phi(t)],
                t=0..30,
                r=0..10,
                p=0..2*Pi,
                phi=0..Pi,
                arrows=cheap,
                dirfield=[5,5,5]
              );

will produce this plot

 

 

 

No...

tomleslie 13428
March 08 2023
2 0

just your misunderstanding of how a plot command "works". When plotting some expression, the plot() command will evaluate the expression at ~200 points (by default), and then produce a "straight line" between adjacent points. The fact that you "see" this as a smooth "curve" is an optical illusion.

In your first example the points are only about 1.8 degrees apart (360/200) by default. In your second example, the points will be 3600000/200 degrees apart (by default) - so 18000degrees or about 50 complete revolutions. - so essentially you end up drawing a straight line from some point on the cicumference of a circle to some other (more or less random) point on the circumference. These lines will fill up the interior of the circle very quickly. You can get around this by increasing the option numpoints drastically and then wait a while for your plot to appear

You appear to get the...

tomleslie 13428
March 02 2023
1 3

@mmcdara 

the styles you want using conventional 'labels' and 'labelfont' options, see the attached

NB the double square brackets around the unit only appear on this site - they do not occur in the original worksheet

> 

  restart:
  with(Units):
  plot( x,
        x=0..1,
        labels=[ typeset(foo*Unit(kg/m^3)),
                 typeset(bar*Unit(m) )
               ],
        labelfont=[times, bold, 14]
      );

Automatically loading the Units[Simple] subpackage
 

  

 
> 

 

> 

 

Download labs.mw

A couple of observations...

tomleslie 13428
March 01 2023
1 3

 

  1. Two of your constraints are on the form s<p<q. Each of these has to be split into two simple inequalities, ie s<p, p<q. This fix is trivial
  2. Both the expression you want to minimise and one of the constraints contain terms such as 91201.82720*(-0.05*i + 0.04)*sqrt(-i), - 912.0182718*i/sqrt(-i). Since 'i' is restricted to range 0..1, the sqrt(-i) factor in both of these terms is going to produce a complex value, which the optimizer cannot handle. ( All optimizers essentially work on the basis of determining whether an expression using the current variable values is greater than or less than the expression with the previous set of variables. However for complex values, the term "greater than or less than" has no meaning. So you have to ask, why do complex values occur?? Is this expected?
  3. If the existence of complex values is unexpected then there must be an error in yur derivation somewhere. These complex valuesappear to arise from the evaluation of the quantities 'Q' and 'Qr'

So many syntax errors...

tomleslie 13428
February 28 2023
0 0

I suggest you read the help page at ?pdsolve/numeric, very carefully

The attached at least executes, although I can't say I'm wildly convinced by the answer which Maple provides

Inline displlay of worksheets on this site appears to be broken - again!

Download pde.mw.

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