Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 163 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Cheat!!!...

tomleslie 13876
January 13 2022
1 1

Since your z-ccordinate is always zero, you don't really have a "3D" curve, you have a 2D curve embedded in a 3D space.

So actually you can  generate/fill the curve by using a 2D parametric plot as shown in the attached

> 

  restart;
  with(plots):
  #Circle of reference
  R0 := 1:
  C0 := t-> < cos(t), sin(t), 0 >:

  #Variable radius
  R := t-> R0*(1 + 1/2*sin(t)):

  #Variable phase
  Dt := t-> Pi/2*sin(2*t):
  
  #CONCAVE Curve
  C := t-> R(t)*C0(t - Dt(t)):
  'C(t)' = C(t):

  #Range
  t1, t2 := 0, 2*Pi:

  #Plot
  GC := spacecurve(C(t),
                   t = t1..t2,
                   color = "Red",
                   thickness = 3,
                   linestyle = solid,
                   scaling = constrained,
                   axes = frame,
                   orientation = [50,40,0]
                 );
  p1:= plot( [ convert(C(t), list)[1..2][],
               t=0..2*Pi
             ],
             scaling=constrained,
             axes=boxed,
             filled=true,
             color=red
           );

 

 
> 

 

> 

 

Download fillcurve.mw

Maybe...

tomleslie 13876
January 09 2022
1 2

you just need to add the 'connect=true' option to the pointplot() command, as in the attached (which for some reason won't display inline here)

malt.mw

In future please use the big green up-arrow in the Mapleprimes toolbar to upload an actual worksheet - it make life so much easier for anyone responding

So...

tomleslie 13876
January 08 2022
0 1

if I strip out all the package loads you don't need, and all the fancy formatting that you don't want/need, and just produce the answer you want - as in the attached - please explain why this isn't sufficient for your purpose

> 

  restart;
  with(VariationalCalculus):
  EulerLagrange(1/2*m*(diff(Y(t), t)^2 + diff(w(Y(t), t), t)^2), t, Y(t));
  convert(%, diff);

{m*((D[1](w))(Y(t), t)*(diff(Y(t), t))+(D[2](w))(Y(t), t))*((D[1, 1](w))(Y(t), t)*(diff(Y(t), t))+(D[1, 2](w))(Y(t), t))-(1/2)*m*(2*(diff(diff(Y(t), t), t))+2*(((D[1, 1](w))(Y(t), t)*(diff(Y(t), t))+(D[1, 2](w))(Y(t), t))*(diff(Y(t), t))+(D[1](w))(Y(t), t)*(diff(diff(Y(t), t), t))+(D[1, 2](w))(Y(t), t)*(diff(Y(t), t))+(D[2, 2](w))(Y(t), t))*(D[1](w))(Y(t), t)+2*((D[1](w))(Y(t), t)*(diff(Y(t), t))+(D[2](w))(Y(t), t))*((D[1, 1](w))(Y(t), t)*(diff(Y(t), t))+(D[1, 2](w))(Y(t), t)))}

  

{m*((eval(diff(w(t1, t), t1), {t1 = Y(t)}))*(diff(Y(t), t))+eval(diff(w(t1, t), t), {t1 = Y(t)}))*((eval(diff(diff(w(t1, t), t1), t1), {t1 = Y(t)}))*(diff(Y(t), t))+eval(diff(diff(w(t1, t), t), t1), {t1 = Y(t)}))-(1/2)*m*(2*(diff(diff(Y(t), t), t))+2*(((eval(diff(diff(w(t1, t), t1), t1), {t1 = Y(t)}))*(diff(Y(t), t))+eval(diff(diff(w(t1, t), t), t1), {t1 = Y(t)}))*(diff(Y(t), t))+(eval(diff(w(t1, t), t1), {t1 = Y(t)}))*(diff(diff(Y(t), t), t))+(eval(diff(diff(w(t1, t), t), t1), {t1 = Y(t)}))*(diff(Y(t), t))+eval(diff(diff(w(t1, t), t), t), {t1 = Y(t)}))*(eval(diff(w(t1, t), t1), {t1 = Y(t)}))+2*((eval(diff(w(t1, t), t1), {t1 = Y(t)}))*(diff(Y(t), t))+eval(diff(w(t1, t), t), {t1 = Y(t)}))*((eval(diff(diff(w(t1, t), t1), t1), {t1 = Y(t)}))*(diff(Y(t), t))+eval(diff(diff(w(t1, t), t), t1), {t1 = Y(t)})))}

 (1)
> 

 

Download EL.mw

Don't understnd the problem...

tomleslie 13876
January 08 2022
1 1

The PhaseSI() command is implemented in the CoolProp() package, and it returns a srting - see the attached.

So what exectly is the problem?

> 

  restart;
  ThermophysicalData:-CoolProp:-PhaseSI(P, 101325, T, 373, Water);
  ThermophysicalData:-CoolProp:-PhaseSI(P, 101325, T, 473, Water);

"liquid"

  

"gas"

 (1)
> 

 

Download CP.mw

Synthesize!...

tomleslie 13876
January 07 2022
1 1

Just pick any three integer values for x, y, z. The three expressions

x*y*z + y*z + y, x*y*z + x*z + z, x*y*z + x*y + x

will then evaluate to integers - whivh provides a possible set of values for a, b, and c

The complicated expression...

tomleslie 13876
January 07 2022
4 0

in your main "for loop" contains the term

t[j - 1]^(gam - 1)

since gam=1/2, this evaluates to

1/sqrt( t[j - 1] )

so when j=1, this evaluates to

1/sqrt(t[0]);

and  t[0] is earlier evaluated to 0, so this  term becomes

1/sqrt(0)

which is undefined - and it is downhill from there!

Well...

tomleslie 13876
January 06 2022
2 2

one can show numerically that it is "plausible". You can play with the value of 's' (>1) in the attached, as well as the number of primes.

 I suspect that provong it analytically would get very difficult - forming the product term for any infinnte series of primes, looks formidable

> 

  restart:
  f:= proc( numprimes, sVal)
           local i,p,prs:=[seq(ithprime(i), i=1..numprimes)]:
           return evalf(eval( mul( 1/(1-1/p^s), p in prs), s=sVal)),
                  evalf(eval( sum( 1/n^s, n=1..infinity), s=sVal))
      end proc:

  f(100, 3);
  f(1000, 3);
  f(10000,3);
  
  f(100, 5);
  f(1000, 5);
  f(10000,5);
           

1.202056602, 1.202056903

  

1.202056902, 1.202056903

  

1.202056903, 1.202056903

  

1.036927755, 1.036927755

  

1.036927755, 1.036927755

  

1.036927755, 1.036927755

 (1)
> 

 

Download euProd.mw

Because...

tomleslie 13876
January 06 2022
1 0

the condition for executing the "while" loop is not satisfied with the values you supply. This is demonstrated in the attached simply by inserting a (temporary) printf() statement. Check the final value printed.

Thus the last statement in your procedure which actually executes is the final "local" declaration, which returns nothing.

Your procedure *may* work for other values of the input arguments.

restart

NULL
f := proc (pp, tm) local TMP1, TMP2, TMP3, TMP4; TMP1 := pp; TMP2 := TMP1-pp^((37-tm)*(5.49*10^(-11)*TMP1^3.88+0.71e-1)/(9.72*10^(-9)*TMP1^3.88+2.3)); TMP3 := 1+5.039*10^(-9)*TMP1^2.88*(TMP1-TMP2)*(37-tm)/(9.72*10^(-9)*TMP1^3.88+2.3)^2; TMP4 := TMP1-TMP2/TMP3; printf("%12.8f %12.8f %12.8f %12.8f %12.8f\n", TMP1, TMP2, TMP3, TMP4, abs(TMP4/TMP1-1)); while abs(TMP4/TMP1-1) <= 0.1e-2 do TMP1 := TMP4; TMP2 := TMP1-pp^((37-tm)*(5.49*TMP1^3.88/10^11+0.71e-1)/(9.72*TMP1^3.88/10^9+2.3)); TMP3 := 1+5.039*TMP1^2.88*(TMP1-TMP2)*(37-tm)/(10^9*(9.72*TMP1^3.88/10^9+2.3)^2); TMP4 := TMP1-TMP2/TMP3; printf("%12.8f %12.8f %12.8f %12.8f %12.8 f\n", TMP1, TMP2, TMP3, TMP4, abs(TMP4/TMP1-1)) end do end proc

f(230, 34)

230.00000000  228.83859720    1.00041038    1.25527390   0.99454229

  

NULL

Download errProc.mw

 

A suggestion...

tomleslie 13876
January 05 2022
1 2

Try inserting the command

Digits:=20;

or even

Digits:=50;

immediately after your 'restart' command and then examine the results - very carefully

Then ask yourself a serious question about just how accurate you would expect two different floating point calculations to be

Use simplfy(expr, zero)...

tomleslie 13876
January 05 2022
3 0

as in the attached

interface(version)

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

 (1)

NULLsqrt(2)*(EllipticK((1/2)*sqrt(2+2*sin(`&phi;__0`)))-EllipticF(sqrt(2)/sqrt(2+2*sin(`&phi;__0`)), (1/2)*sqrt(2+2*sin(`&phi;__0`))))

2^(1/2)*(EllipticK((1/2)*(2+2*sin(phi__0))^(1/2))-EllipticF(2^(1/2)/(2+2*sin(phi__0))^(1/2), (1/2)*(2+2*sin(phi__0))^(1/2)))

 (2)

2*EllipticF(sqrt(sin(`&phi;__0`))/sqrt(sin(`&phi;__0`)+1), I*sqrt(-sin(`&phi;__0`)^2+1)/(sin(`&phi;__0`)-1))/sqrt(1-sin(`&phi;__0`))

2*EllipticF(sin(phi__0)^(1/2)/(1+sin(phi__0))^(1/2), I*(-sin(phi__0)^2+1)^(1/2)/(sin(phi__0)-1))/(1-sin(phi__0))^(1/2)

 (3)

simplify(eval(2^(1/2)*(EllipticK((1/2)*(2+2*sin(phi__0))^(1/2))-EllipticF(2^(1/2)/(2+2*sin(phi__0))^(1/2), (1/2)*(2+2*sin(phi__0))^(1/2))), `&phi;__0` = .1), zero)

.448035371*2^(1/2)

 (4)

simplify(eval(2*EllipticF(sin(phi__0)^(1/2)/(1+sin(phi__0))^(1/2), I*(-sin(phi__0)^2+1)^(1/2)/(sin(phi__0)-1))/(1-sin(phi__0))^(1/2), `&phi;__0` = .1), zero)

.6336176952

 (5)

op(.6336176952)

6336176952, -10

 (6)

whattype(6336176952, -10)

exprseq

 (7)

0.*I;

0.*I

 (8)

whattype(0.*I); simplify(0.*I, zero); whattype(%)

float

 (9)

0.*I

0.*I

 (10)

whattype(0.*I)

complex(extended_numeric)

 (11)

NULL

Download zeroI.mw

Consider the zip() command...

tomleslie 13876
January 04 2022
1 2

There is nothing intrinsically *wrong* with your approach of using elementwise operations, but (as Acer has pointed out), this can become difficult to read/debug. In Maple there is often more than one way to do it. For this particular problem, I would probably use the zip() command as in the attached.

> 

  V1 := Vector[column](5, [31.6, 30.3, 32.6, 35, 32.7]):
  V2 := Vector[column](5, [225.9, 161.7, 147.2, 215, 207]):
  f1:= (x,y)-> x*10^(0.032077*(y-37));
  f2:= (x,y)-> x*10^((5.49*10^(-11)*x^3.88+0.071)/(9.7*10^(-9)*x^3.88+2.30)*(y-37));
  zip( f1, V2, V1);
  zip( f2, V2, V1);

f1 := proc (x, y) options operator, arrow; x*10^(0.32077e-1*(y-37)) end proc

  

f2 := proc (x, y) options operator, arrow; x*10^((5.49*(1/100000000000)*x^3.88+0.71e-1)*(y-37)/(9.7*(1/1000000000)*x^3.88+2.30)) end proc

  

Vector(5, {(1) = 151.6004292, (2) = 98.58120443, (3) = 106.3577179, (4) = 185.4746072, (5) = 150.6743218})

  

Vector[column](%id = 36893488148091976628)

 (1)
> 

 


 

Download zipV1V2.mw

 

Your code...

tomleslie 13876
December 28 2021
3 6

and question confuse me.

You state that Nb must be equal to Nt but in your code these are varied independently - why?

You state that you want to plot -diff(theta(eta),eta) vs Bi, but in your code, Bi is fixed at 0.1 - why?

The diagram you show, plots -diff(theta(eta),eta) evaluated at theta=0 vs Bi, is this really what you want? If it is then the attached shows how to obtain and plot -diff(theta(eta),eta) evaluated at theta=0 vs Bi for Nb=Nt=0.2 and Bi=0..5

> 

restart
> 

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 = 10, Le = 10, 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:
  plot
  ( getRes( 0.2, x),
    x=0..5
  );
               

 
> 

 

> 

 

Download odeProb.mw

Hmmm...

tomleslie 13876
December 24 2021
2 1

Let us consider your original equation, pasted from your original post

restart;
eq1 := diff(u(u, t), t) = A[0] + A[1]*cos*omega*t + Beta[1]*[diff(u(u, r), r $ 2) + diff(u(u, r), r)/r];

Questions

  1. Exactly what is the dependent function - u(u, t) or u(u, r) - because you can't have both.
  2. Do not use the same name as both a dependent function and an inddependent variable - as in u(u, t)
  3.  Maybe the dependent function is actually u(r, t)?
  4. When you write cos*omega*t, did you actually mean cos(omega*t) ?
  5. Do not use square brackets ( is [] ) for simple term grouping such as [diff(u(u, r), r $ 2) + diff(u(u, r), r)/r]. Square and curly (ie {}) brackets have special meanings in Maple and only round brackets (ie '()' ) should be used for term grouping

It wuld be helpful if you posted a worksheet - using the big green up-arrow in the Mapleprimes toolbar

My $0.02...

tomleslie 13876
December 23 2021
0 0

using your data points - although the way you have uploaded the prtoblem I only have access to the first 10 points

The following will produce the plot you want (for the first 10 points), and (in principle) will do so even if you supply a 5681x3 matrix - although this will generate 5680 frames - so you might have to wait a while. I agree with mmcdara that a better option might be to plot more than one "new" point per frame. Adjusting the number of frames on the basis of the divisors of 5681 would be pretty trivial.

> 

  restart;
  with(plots):

  pts:= Matrix
        ( [ [-0.604017052608010,-0.264977766339890,0.751869516706500],
            [-0.613978487863800,-0.249414876060040,0.749375367945900],
            [-0.623323756089740,-0.233553456286040,0.747039833383040],
            [-0.632047442843230,-0.217415623800280,0.744865677145030],
            [-0.640144714635430,-0.201023189254840,0.742855445252480],
            [-0.647611273173430,-0.184397665673180,0.741011468595750],
            [-0.654443312425730,-0.167560280356470,0.739335865792190],
            [-0.660637478673220,-0.150531989900910,0.737830545893180],
            [-0.666190833687370,-0.133333498035410,0.736497210912830],
            [-0.671100821158590,-0.115985275993310,0.735337358152890]
          ]
        ):
  animate
  ( pointplot3d,
    [ 'pts[i-1..i, 1..3]',
      color=blue,
      style=pointline,
      symbol=solidsphere,
      symbolsize=16,
      thickness=4
    ],
    i=[$ 2..op([1,1], pts)],
    trace = op([1,1], pts)
  );

 
> 

 

NULL

Download amin3.mw

the problem...

tomleslie 13876
December 22 2021
2 1

is that 0*Unit(anyUnit) automatically simplifies to 0 (without a unit) whihc occasionally causes problems.

One workaround for your specific case is shown in the attached

> 

  restart;
  with(Units[Simple]):
  with(plots):
  diffdiameter := (conc, number) -> -number*diameter*exp(-conc/Unit(mol/kg));
   
  diameter := 2*Unit(m);

  evalf(diffdiameter(1*Unit(mol/kg), 6));
  evalf(diffdiameter(1*Unit(mol/kg), 0));

  convert(evalf(diffdiameter(Unit(mol/kg), 6)), unit_free);
  convert(evalf(diffdiameter(Unit(mol/kg), 0)), unit_free);

  display
  ( [ plot
      (  convert(evalf(diffdiameter(conc, 0)), unit_free),
         conc = 0 .. 4,
         color=blue
      ),
      plot
      ( convert(evalf(diffdiameter(conc, 6)), unit_free),
        conc = 0 .. 4*Unit(mol/kg),
        color=red
      )
    ],
    title = "derivative Diameter with conc",
    labels = ["concentration / mol/kg", "deriv diameter / m"],
    color = ["blue", "red"],
    legend = ["bare diameter", "PLUS diameter"],
    labeldirections = ["horizontal", "vertical"],
    titlefont = [Helvetica, bold, 16],
    axesfont = [Helvetica, 14],
    labelfont = [Helvetica, 14],
    axes = boxed
 );

proc (conc, number) options operator, arrow; Units:-Simple:-`-`(Units:-Simple:-`*`(number, diameter, Units:-Simple:-exp(Units:-Simple:-`-`(Units:-Simple:-`*`(conc, Units:-Simple:-`/`(Unit(Units:-Simple:-`*`(mol, Units:-Simple:-`/`(kg))))))))) end proc

  

2*Units:-Unit(m)

  

-4.414553294*Units:-Unit(m)

  

0.

  

-4.414553294

  

0.

  

 
> 

 

Download unitProb.mw

First 34 35 36 37 38 39 40 Last Page 36 of 207