mmcdara

7891 Reputation

22 Badges

9 years, 51 days

MaplePrimes Activity


These are answers submitted by mmcdara


Orignal reply deleted by the author.

Revised version: Updated.mw

 


From the reference you give:

 

restart

a := 1;

h := 1;

# pA := [0, 0, 0]; # useless

pB := [a, 0, 0];

pC := [a, a, 0];

pD := [0, a, 0];

pS := [0, 0, h];

1

 

1

 

[1, 0, 0]

 

[1, 1, 0]

 

[0, 1, 0]

 

[0, 0, 1]

(1)

local D:
symbols := [B, C, D, S]:

use geom3d in
  seq(point(s, p||s), s in symbols);
  plane(SBC, [S, B, C]);
  plane(SCD, [S, C, D]);
  sbc    := < NormalVector(SBC) >;
  scd    := < NormalVector(SCD) >;

  # from the reference you give :
  # "the dihedral angle, φ... satisfies 0 ≤ φ ≤ π/2"
  # So 2*π/3 and 3*π/4 cannot be solutions as you say.
  
  cosphi := abs(Statistics:-Correlation(sbc, scd));
end use:

phi := identify(arccos(cosphi))

(1/3)*Pi

(2)

 


 

Download Dihedral_Angle.mw


As you've already noticed, there's nothing immediate..
The trick to derive equations 9, 11 and 13 is to multiply both sides of equation 7 by G' and integrate wrt eta (without without forgetting the "constant" K)

(Done with Maple 2015)
 

restart

eq_7 := diff(G(eta), eta$2) + sigma*G(eta) = nu

diff(diff(G(eta), eta), eta)+sigma*G(eta) = nu

(1)

 

sigma < 0

 

negsol := rhs( ( dsolve(eq_7) assuming sigma < 0) ):
negsol := expand(convert(negsol, trig)):


CS      := [C, S]:
replace := convert(indets(negsol, function), list) =~ CS:
negsol  := collect(eval(negsol, replace), CS):
B1B2    := solve({coeff(negsol, C)=B1, coeff(negsol, S)=B2}, [_C1, _C2]):

eq_8 := eval(eval(negsol, %[]), (rhs=lhs)~(replace));

B1*cosh((-sigma)^(1/2)*eta)+B2*sinh((-sigma)^(1/2)*eta)+nu/sigma

(2)

use IntegrationTools in

  # multiply both member of eq_7 by G'

  expand~(eq_7*~diff(G(eta), eta));
 
  # Integrate over eta

  map((Expand@Int), %, eta);
 
  # evaluate

  value~(%);
 
  # divide each term by G^2 (assuming G <> 0)

  expand(% /~ G(eta)^2);
 
  # isolate (G'/G)^2

  ref := isolate(%, (diff(G(eta), eta) / G(eta))^2);
 
  # add an integration "constant"

  ref := lhs(ref) = rhs(ref) + __K;
 
  # plug in ref the expression of G(eta) given by eq_8

  aux := eval(ref, G = (eta -> eq_8));

  # solve for __K

  K := simplify(solve(aux, __K));

  # Identify the denominator of K

  den := simplify(denom(K), {eq_8=G(eta)});

  # Rewrite K according to theprevious result

  K := numer(K)/den;

  # Finally

  eq_9 := eval(ref, __K=K);

end use:

eq_9;

(diff(G(eta), eta))^2/G(eta)^2 = 2*nu/G(eta)-sigma+(B1^2*sigma^2-B2^2*sigma^2-nu^2)/(sigma*G(eta)^2)

(3)

 

sigma > 0

 

eq_10 := eval( rhs( ( dsolve(eq_7) assuming sigma > 0) ), [_C1=B1, _C2=B2])

sin(sigma^(1/2)*eta)*B2+cos(sigma^(1/2)*eta)*B1+nu/sigma

(4)

use IntegrationTools in
 
  # plug in ref the expression of G(eta) given by eq_8

  aux := eval(ref, G = (eta -> eq_10));

  # solve for __K

  K := simplify(solve(aux, __K), size):
  K := simplify(numer(K))/denom(K);

  # Identify the denominator of K

  den := simplify(denom(K), {eq_10=G(eta)});

  # Rewrite K according to theprevious result

  K := numer(K)/den;

  # Finally

  eq_11 := eval(ref, __K=K);

end use:

eq_11

(diff(G(eta), eta))^2/G(eta)^2 = 2*nu/G(eta)-sigma+(B1^2*sigma^2+B2^2*sigma^2-nu^2)/(sigma*G(eta)^2)

(5)

 

sigma = 0

 

zerosol := dsolve(eval(eq_7, sigma = 0)):
eq_12   := rhs(eval(zerosol, {_C1=B1, _C2=B2}));

(1/2)*nu*eta^2+B1*eta+B2

(6)

use IntegrationTools in

  # multiply both member of eq_7 by G'

  expand~(eval(eq_7, sigma=0)*~diff(G(eta), eta));
 
  # Integrate over eta

  map((Expand@Int), %, eta);
 
  # evaluate

  value~(%);
 
  # divide each term by G^2 (assuming G <> 0)

  expand(% /~ G(eta)^2);
 
  # isolate (G'/G)^2

  ref := isolate(%, (diff(G(eta), eta) / G(eta))^2);
 
  # add an integration "constant"

  ref := lhs(ref) = rhs(ref) + __K;
 
  # plug in ref the expression of G(eta) given by eq_8

  aux := eval(ref, G = (eta -> eq_12));

  # solve for __K

  K := simplify(solve(aux, __K), size):
  K := simplify(numer(K))/denom(K);

  # Identify the denominator of K

  den := simplify(denom(K), {eq_12=G(eta)});

  # Rewrite K according to theprevious result

  K := simplify(numer(K)/den);

  # Finally

  eq_13 := eval(ref, __K=K);

end use:

eq_13

(diff(G(eta), eta))^2/G(eta)^2 = 2*nu/G(eta)+(B1^2-2*B2*nu)/G(eta)^2

(7)

 


 

Download This_way.mw



It's always annoying when the question you ask is not answered nor commented. 
Does this mean that nobody has any answers, or that the question isn't even worth spending time answering because it's completely irrelevant?
I do not consider your question lacks interest.

I'm afraid there is no solution to your problem, even in the far simpler 2D case.

Here are a few ideas ( a white elephant?) to rotate a text
Aligned.mw


This is bery basic stuff that must be enhanced to give acceptable results: increase the size of the text in t:=... (for instance font=[Helvetica, 50] for a better definition), adjust the threhold in function Tf , adjust the gridsize in densityplot, ...).

Hope this will help

restart

# Access to the help page

?sum

# from ,sum help page, example 4 (Maple 2015)

sum(1/k!,k=0..infinity);
evalf(%);

exp(1)

 

2.718281828

(1)

# A first example invoking the exponential function

sum(exp(-k), k=0..infinity);

-1/(exp(-1)-1)

(2)

# Use of the Inert form ('Sum'):

Sum(exp(-k), k=0..infinity) = sum(exp(-k), k=0..infinity);

Sum(exp(-k), k = 0 .. infinity) = -1/(exp(-1)-1)

(3)

# Be careful: functions 'add' and 'sum' do mean the same thing.
# This is a common error.
#
# To understand the difference between them: 'add' stands for 'addition' while
# 'sum' stands for... 'sum'
# (see these terms in Maple's dictionary: within the "Maple xxx Help" window: in the
# left pane click on the drop-down box labeled "Page Types",  unclick "Select all"
# and select only "Dictionary").
#
# Once done type addition in the top-left field and "Search": the right page should
# the definition of the term addition and refer to the function contain 'add'.
#
# Next type sum instead of addition and read carefuly the second item in the right pane:
#    2.  The limit of the sequence  of partial sums of the first n terms of an infinite series ,
#        as n tends to infinity.
#

 

 

Download sum.mw

Can you please keep this as a comment and avoid converting it into an answer, there is no reason to do so.

After your dsolve command type this to see that sys still contains an unknown quantity named e for which no ic/bc condition is given, so the error you get:

indets(sys1, name) minus {eta}
                              {e}

Cure: fix the value of e before executing dsolve or do

dsolve(eval(sys, e=some_numeric_value), numeric)

For instance

sol1 := dsolve(eval(sys1, e = 1), numeric);
           sol1 := proc(x_bvp)  ...  end;

the answer is NO.
But you can easily use Eclipse, GitHub or any other versioning tool provided you use mpl files and not mw files.

If you have several mw files it will take some time to convert them all by hand (open the mw file ans export it as mpl) and I don't know if these actions can be done programmatically.
Once the mpl files are created I advice you to use some tool which enables syntaxic coloring.
For my part, I use GitHub with Notepad++ but I know a coleague of mine uses Eclipse.

A hint: search "versioning" on this site, if I'm not mistaken I asked a similar question years ago (maybe from my professional account sand15).

Can you please keep this as a comment and avoid converting it into an answer, there is no reason to do so.

restart:

with(geometry):

with(plots):

 

Bl := color = black;

y0 := x -> -ln(1 - exp(-x)):

y1 := x -> -ln(-1 + exp(-x)):

y2 := x -> -ln(1 + exp(-x)):

 

p  := plot(y0(x), x = 0.02 .. 4, scaling = constrained, color = blue):

p1 := plot(y2(x), x = -4 .. 4, scaling = constrained, color = green):

p2 := plot(y1(x), x = -4 .. 0, scaling = constrained, color = red):

display({p, p1, p2}, view = [-4 .. 4, -4 .. 5]);

color = black

 

 

# First way

plots:-display(
  plots:-inequal({y >= y2(x), y <= y1(x)}, x=-4..0, y=-4..4),
  plots:-inequal({y >= y2(x), y <= y0(x)}, x=0..4, y=-4..4)
);

 

# A less efficient way but having the characteristic function of the blue domain
# might be interesting.

interior := (x, y) -> piecewise(
                        x >=0
                        , piecewise(y < y0(x) and y > y2(x), 1, 0)
                        , piecewise(y < y1(x) and y > y2(x), 1, 0)
                      ):

plots:-contourplot(
  interior(x, y)
  , x=-4..4, y=-4..4
  , contours=1
  , grid=[200, 200]    # use grid=[400, 400] for a better rendering
  , filled=true
)

 

Area := int(int(1, y=y2(x)..y0(x)), x=0..4) + int(int(1, y=y2(x)..y1(x)), x=-4..0);
evalf(Area);

(1/6)*Pi^2-dilog(1-exp(-4))+dilog(1+exp(-4))-dilog(exp(4))-4*ln(-1+exp(4))-dilog(1+exp(4))

 

4.861536918

(1)
 

 

Download Amoeba.mw

A few workarounds colors.mw
 

restart


Trigonometric functions are defined over radians, not degrees nor grades.

An elementary function cannot contain a dimension, unless it is a trigonometric function and
the dimension is 'radian'.

# Example 1: radians times meters

3*Units:-Unit('radians') * 2*Units:-Unit('meters');
simplify(%)

6*Units:-Unit('rad')*Units:-Unit('m')

 

6*Units:-Unit(('m')^2/('m')('radius'))

(1)

# Example 2: The 'ln' case

ln(3.14);

1.144222800

(2)

ln(3.14*Units:-Unit('m'));
simplify(%);

ln(3.14*Units:-Unit('m'))

 

Error, (in Units:-Standard:-ln) the function `Units:-Standard:-ln` cannot handle the unit m with dimension length

 

ln(3.14*Units:-Unit('radians')),
simplify(%)

Error, (in Units:-Standard:-ln) the function `Units:-Standard:-ln` cannot handle the unit m with dimension length

 

# Example 3: The trigonometric case
#
# Here the argument of 'tan' has 'radians' unit

tan(Pi/4);
tan(convert(45, 'units', 'degrees', 'radians'));

# here it has 'degrees' unit

tan(45*Units:-Unit('degrees'));

# note that simplify returns the correct result, likely because simplify
# converts back the degrees into radians (?)
simplify(%)

1

 

1

 

tan(45*Units:-Unit('arcdeg'))

 

1

(3)

# Your example

eq := x -> -log(tan(x/2)-3)/5 + ln(3*tan(x/2)+1)/5

proc (x) options operator, arrow; -(1/5)*log(tan((1/2)*x)-3)+(1/5)*ln(3*tan((1/2)*x)+1) end proc

(4)

# First way

xdeg := 10:
xrad := convert(xdeg, 'units', 'degrees', 'radians');

eq(xrad);
evalf(%);

(1/18)*Pi

 

-(1/5)*ln(tan((1/36)*Pi)-3)+(1/5)*ln(3*tan((1/36)*Pi)+1)

 

-.1671897534-.6283185308*I

(5)

# Second way

xdeg := 10*Units:-Unit('degrees'):
eq(xdeg);
simplify(%);
evalf(%);

-(1/5)*ln(tan(5*Units:-Unit('arcdeg'))-3)+(1/5)*ln(3*tan(5*Units:-Unit('arcdeg'))+1)

 

-(1/5)*ln(tan((1/36)*Pi)-3)+(1/5)*ln(3*tan((1/36)*Pi)+1)

 

-.1671897534-.6283185308*I

(6)
 

 

Download degrees.mw


which works for polynomial expressions of whose all functions depend on the same indeterminate.

For instance

expr = x(t)*diff(z(t), t)^2,
expr = diff(x(t), t$2)+x(t)*diff(z(t), t$3),
...

but not

expr = x(t)*z(u)^2,
expr = x(t)*diff(cos(u(t)), t),
...

 
 

restart

ddf := proc(F, {n::posint:=1})

  local f, df, idd, idf, rules, i, deg, opi:
  f  := eval(F, indets(F, name) =~ indets(F, name)(__t));
  df := diff(f, __t$n);

  idd, idf := selectremove(has, indets(df), diff);
  idd := convert(idd, list);
  idf := convert(idf minus {__t}, list);

  rules := NULL:
  for i in idd do
    if is(op([0, 0, 0], convert(i, D)) = `@@`) then
      deg := op(2, op([0, 0], convert(i, D)));
      opi := op([0, 1], convert(i, D));
    else
      deg := 1:
      opi := op([1, 0], i);
    end if:
    rules := rules, i = cat(`#mrow(mi("`, opi, `")`, seq(`,mi("&#x27;")`, j=1..deg), `)`);
  end do:
  for i in idf do
    rules := rules, i = cat(`#mi("`, op(0, i), `")`);
  end do;
  
  eval(df, [rules])
end proc:

 

ddf(z^2+y);
ddf(z^2+y, n=1);
ddf(z^2+y, n=2);
ddf(t, n=4);

2*`#mi("z")`*`#mrow(mi("z"),mi("&#x27;"))`+`#mrow(mi("y"),mi("&#x27;"))`

 

2*`#mi("z")`*`#mrow(mi("z"),mi("&#x27;"))`+`#mrow(mi("y"),mi("&#x27;"))`

 

2*`#mi("z")`*`#mrow(mi("z"),mi("&#x27;"),mi("&#x27;"))`+2*`#mrow(mi("z"),mi("&#x27;"))`^2+`#mrow(mi("y"),mi("&#x27;"),mi("&#x27;"))`

 

`#mrow(mi("t"),mi("&#x27;"),mi("&#x27;"),mi("&#x27;"),mi("&#x27;"))`

(1)

# Your result: y*z*(2*z*z'+y')*y''+6*y^2*x''''

expr := y * z * ddf(z^2+y, n=1) * ddf(y, n=2) + 6 * ddf(x, n=4) *y^2;

y*z*(2*`#mi("z")`*`#mrow(mi("z"),mi("&#x27;"))`+`#mrow(mi("y"),mi("&#x27;"))`)*`#mrow(mi("y"),mi("&#x27;"),mi("&#x27;"))`+6*`#mrow(mi("x"),mi("&#x27;"),mi("&#x27;"),mi("&#x27;"),mi("&#x27;"))`*y^2

(2)

 


 

Download An_alternative.mw

 



Could this help you Fireworks?

Another exemple of trajectory fading is provided in Jonglage_4.mw.
Here is an animation of a ball bouncing within a triangular shape: at the current time T the trajectory is plotted between T-dt and T, a colorscheme whose first color is equal to the background color mimics a fadding trajectory.

Based on Thomas' example  you would obaiin this


Thomas_with_fading.mw

 

It is not because the display of den1 contains terms involving t that den1 INDEED depends on t.
This is a "side effect" of simplify.


Functions_of_t := print~(select(has, indets(simplify(eval(den1, para), size), function), t)):
                        exp(-I exp(I t))

                            exp(I t)

                        exp(I exp(I t))

               /1   /   (1/2)  (1/2)             \\
            exp|- I \x 3      4      - 4 exp(I t)/|
               \4                                 /

               /1   /   (1/2)  (1/2)             \\
            exp|- I \x 3      4      + 4 exp(I t)/|
               \4                                 /


 Functions_of_t := select(has, indets(expand(eval(den1, para)), function), t);
                               {}

numericden_mmcdara.mw

It is easy to verify that G(x, y) = G(y, x).
Then the minimum is to be search along the line y=x.

# Verify that G(x, y)
simplify(G(x, y)-G(y, x))
                               0

# Define a univariate quantity T
T := expand(G(x, x));
s := [fsolve(diff(T, x))];
                  s := [0., .8073218015, 1.015827227]

# Which one is a maximizer/minimizer ?
 eval~(diff(T, x$2), x=~s);
                 [-5.600000000, 30.91991289, -88.23416528]

# then s[2] is the minimizer and the minimum is equal tools

G(s[2], s[2])
                          -0.397001070

# With minimize
minimize(T, x=0..1, location=true);
      -0.3970010694, {[{x = 0.8073218015}, -0.3970010694]}

(In Statistics a model is said to be linear if it writes as a linear combination of its parameters, not of its regressors [here lambda])

Here obj simply writes C__10*U + C__10*V where both U and V are known functions of lambda.
Then your problem is linear wrt C__10 and C__10LINEAR.mw 


Linearization is by far the simplest way to solve the problem (you can even find the solution by hand).

First 8 9 10 11 12 13 14 Last Page 10 of 65