mmcdara

5274 Reputation

17 Badges

7 years, 117 days

MaplePrimes Activity


These are answers submitted by mmcdara

 

One way which introduces only a slight modification to your code is to replace the line

eq1 := Qbar = InvT . Q . R . T . InvR

by the red line below

assign([entries(Qbar, nolist)] =~ [entries(InvT . Q . R . T . InvR, nolist)]);

plot(Qb11, p = 0 .. 9): # fails for Qb11 depends on 3 variables:

Warning, expecting only range variable p in expression (.2339357430e12*cos(p)^2+.2319277108e11*s^2)*cos(p)^2+(.2319277108e11*cos(p)^2+s^2*Q22)*s^2+.2868e11*s^2*cos(p)^2 to be plotted but found names [Q22, s]

# check this
indets(Qb11, name);

                          {Q22, p, s}

# do this, for instance, to evaluate Qb11 for siome values of Q22 and s
plot(eval(Qb11, [Q22=1, s=1]), p=0..9)

 

with recent Maple versions

numer(a) %/ denom(a);


For older versions (2015 for instance)
 

with(InertForm):
n := convert(numer(a), string):
d := convert(denom(a), string):
Parse( cat("(", n, ")/", d) ):
Display(%);

 


Hints for Problem 7
An equivalence relation E has 3 properties R, S and T (which are?).

  • one of them (R) is obvious
  • one (S) comes from the commutativity of the addition
  • the last (T) involes 3 couples C, C' and C"
    • write C E C' and C' E C""
    • try to combine the two equalities to obtain C E C" (needs a simplification and uses properties of the addition)
       
  • In fact, let F a relation which has the same properties of the addition, then  (a, b) E (c, d) defined by  (a+d) F (b+c) is an equivalence relation.

 

For the three other problems look at some hints here
pbs.pdf


It's clear that D is at the intersection of the bisector of angle A and the line which passes through B and C.
Thus its coordinates and a parametric representation of AD

coords__D := ([-6, 1,-1] +~ [2,-3, 7]) /~ 2;
line__AD  := [x, y, z] =~ [-2, 3, 5]*~s +~ coords__D*~(1-s);

 

But as I'm fond of the geom3d package, here is a step by step wauy yo use it (uselesss for this problem but it could help you for more complex ones). 

restart:
with(geom3d):

point(A, -2, 3, 5):
point(B, -6, 1,-1):
point(C,  2,-3, 7):

line(AB, [A, B]):
line(AC, [A, C]):

# parametric equations of AB and AC

eq__AB := [x, y, z] =~ Equation(AB, s):
eq__AC := [x, y, z] =~ Equation(AC, s):

# AD is the bisector of angle BAC:
# parametric equation of AD

eq__AD := (eq__AB +~ eq__AC) /~ 2;
               [x = -2, y = 3 - 4 s, z = 5 - 2 s]

# implicit equation of AD 
map(u -> if depends(u, s) then solve(u, s) else rhs(u) end if, eq__AD):
eq__AD__imp := add(%);
                         5   1     1  
                         - - - y - - z
                         4   4     2  

# coordinates of D
local D:
line(BC, [B, C]):
line(AD, rhs~(eq__AD), s):
intersection(D, BC, AD):
detail(D)
          Geom3dDetail(["name of the object", D], 
            ["form of the object", point3d], 
            ["coordinates of the point", [-2, -1, 3]])

# visualize
plots:-display(
  PLOT3D(
    POINTS([coordinates(A)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(B)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(C)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(D)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 0, 0, 1)),
    TEXT(1.1*~coordinates(A), "A", FONT(Helvetica, bold, 14)),
    TEXT(1.1*~coordinates(B), "B", FONT(Helvetica, bold, 14)),
    TEXT(1.1*~coordinates(C), "C", FONT(Helvetica, bold, 14)),
    TEXT(1.2*~coordinates(D), "D", FONT(Helvetica, bold, 14), COLOR(RGB, 0, 0, 1)),
    POLYGONS([coordinates(A), coordinates(B), coordinates(C)], COLOR(RGB, 0.9$3))
  ),
  plot3d([rhs~(eq__AD)[]], s=0..1, color=blue, thickness=5)
)



with_geom3d.mw

I'm not going to risk a stiff neck to read this but it sounds like a fixed point method.
Here is an example, adapt it to your own function f


 

restart:

with(plots):

# example

f := x -> surd(x, 3):
g := x -> surd(x, 4) / 2:

pf := plot(f(x), x=0..1.3, color=blue , thickness=2, gridlines=true):
pg := plot(g(x), x=0..1.3, color=green, thickness=2, gridlines=true):
pl := plot( x  , x=0..1.3, color=black):

start := eval(x, solve({g(x)=x, x > 0}));

(1/4)*2^(2/3)

(1)

tu := `#mo("▲")`;
tr := `#mo("►")`;

`#mo("▲")`

 

`#mo("►")`

(2)

iter := proc(p, eps)
  local q := copy(p):
  local L, r:
  L := q, [q[1], f(q[1])];
  r := [f(q[1]), f(q[1])];
  L := L, r;
  while add(q-~r)^~2 > eps^2 do
    q := r:
    L := L, [q[1], f(q[1])];
    r := [f(q[1]), f(q[1])];
    L := L, r;
  end do:
  return L
end proc:

path := [iter([evalf(start), 0], 0.05)]:
display(
  pf, pg, pl,
  plot(path, linestyle=3),
  seq(textplot([op((path[k]+~path[k+1])/~2), tr]), k in [seq(2..numelems(path), 2)]),
  textplot([op(([start$2]+~path[2])/~2), tu]),
  seq(textplot([op((path[k]+~path[k+1])/~2), tu]), k in [seq(3..numelems(path)-1, 2)]),
  size=[800, 800]
)

 

 


 

Download fixed_point.mw

Look to the first line in the attached fie, this could also interest you

 

For lazzy person only

(PS: you will probably find everything you need in the Finance package; unfortunately, English not being my mother tongue and me not being at all familiar with financial terms, I used a summation to get the result. I have no doubt that if Finance has no secret for you you will find in this package the procedure that answers your question)

restart:
account_earns := 7.75/100;
years         := 30;
investment    := 4000;
                         0.07750000000
                               30
                              4000
add(Finance:-futurevalue(investment, account_earns, years-k), k=0..years-1)
                                      5
                        4.664152587 10 

But you must have been taught how to calculate this?
Does the term Geometric_series ring any bell?
 

Answer to your first question

This is a simple proposal (I used Maple 2015, newer versions produces more beautiful graphs).
I did not plot a ball at the vertices, but draw the path (sequence of edges) from the root of the tree to a termibal leaf).
A ball can be drawn (as you requested) but we need do know if you just want it to be displayed when it arrives on a vertex (simple) or if you want also bo follow it along the edges (more complex).
 

restart

with(GraphTheory):
with(RandomGraphs):
with(plots):

# an example of a random tree graph

T := RandomTree(20):
DrawGraph(T)

 


Part 1: step by step explanation

# leaves (vertices of degree 1)

Leaves := map(v -> if Degree(T, v)=1 then v end if, Vertices(T)[2..-1])

[2, 4, 6, 7, 10, 11, 15, 16, 18, 19, 20]

(1)

# departure vertex

U := 1:

# random selection of the arrival leave

V := combinat:-randperm(Leaves)[1];

# shortest path from departure to arrival vertices

P := ShortestPath(T, U, V);

# coordinates of the vertices in P

X := GetVertexPositions(T)[P]

4

 

[1, 4]

 

[[.5000000000, 1.], [0., .8000000000]]

(2)

# parametric description of the trajectory

traj := CurveFitting:-BSplineCurve(X, t, order=1):

# The last operator (t=0..something) is to be corrected

traj := subsop(3=(t=0..numelems(P)-1), traj):

# Tree + animated path
# example

dg := DrawGraph(T):
animatecurve(traj, background=dg, color=red, thickness=3, frames=numelems(P)):


Part 2: all of this within a procedure

f := proc(T::Graph)
  local U, Leaves, V, P, X, dg, traj:
  U      := 1:
  Leaves := map(v -> if Degree(T, v)=1 then v end if, Vertices(T)[2..-1]);
  V      := combinat:-randperm(Leaves)[1];
  P      := ShortestPath(T, U, V());
  X      := GetVertexPositions(T)[P];
  traj   := CurveFitting:-BSplineCurve(X, t, order=1):
  traj   := subsop(3=(t=0..numelems(P)-1), traj):
  dg     := DrawGraph(T):
  animatecurve(traj, background=dg, color=red, thickness=3, frames=numelems(P));
end proc:

f(T)

 

 

 

Download Animated_path_along_a_tree.mw


Maybe you could be interested in plotting different random paths?
Something like this....

f := proc()
  local U, Leaves, V, P, X, T2, dg, traj:
  U      := 1:
  Leaves := map(v -> if Degree(T, v)=1 then v end if, Vertices(T)[2..-1]);
  V      := combinat:-randperm(Leaves)[1];
  P      := ShortestPath(T, U, V());
  X      := GetVertexPositions(T)[P];
  traj   := CurveFitting:-BSplineCurve(X, t, order=1):
  traj   := subsop(3=(t=0..numelems(P)-1), traj):
  T2     := CopyGraph(T):
  HighlightVertex(T2, P, red); # highlight the vertices the ball hits
  dg     := DrawGraph(T2):
  display(dg, plot(traj, color=red, thickness=3, frames=numelems(P)));
end proc:

N := 10:
#animate(f, [], n=1..N, frames=N)

 

Once you know that your serie is the Taylor expansion of  1/(1+exp(x-5*t))^2 at t=0 (see here 232686-INFINITE-SERIES-SUM)
the answer is obvious

plot3d(1/(1+exp(x-5*t))^2, x=-10..10, t=-2..2, labels="x", "t", "u")

 


It's obvious that your serie can be written as a geometric serie (negative if lambda < 0, alternated if lambda > 0)

restart:
a := lambda*exp(t):
q := (-1)*lambda*(exp(t)-1):
U := (t, n) -> add(a*q^k, k=0..n):

#check
U(t, 4);

# let's set define S(t, lambda, n) this way
S := (t, lambda, n) -> lambda*exp(t) * add(((-1)*lambda*(exp(t)-1))^k, k=0..n):

# S(t, lambda, t, n) converges iif abs(q) < 1
# assuming lambda > 0 S converges iif abs(q) < 1 that is if 0 < lambda < 1/(exp(t)-1)
# Examples

T := 2.0:
plot([seq([k, S(T, 0.95/(exp(T)-1), k)], k=1..200)]);
plot([seq([k, S(T, 1.05/(exp(T)-1), k)], k=1..200)]);

So, what is lambda?

If the serie converges its limit is 

 a * 1/(1-q);
                         lambda exp(t)     
                    -----------------------
                    1 + lambda (exp(t) - 1)

This can be obtained by elementary calculus: the sum of a (converging)  geometric serie is  limit(a*(1-q^n)/(1-q) , n=+oo) = a/(1-q) for, because of the convergence condition abs(q)<1, q^n vanishes to 0.

after the line sys := ... write

indets(sys, name);
                        {a, b, c, d, x}

You have 4 equations and you claime having 3 unknowns.
Thus you miss one (a?) ... or maybe are you looking for another type of solution (the one that would minimize some functionnal of the 4 equations?

I took the liberty to consider a is an unknown (thus x is the only input).
A numeric solution can be obtained that way.

sys := unapply(sys, x):
# example for x=1
fsolve(sys(1))
    {a = -0.2722166549, b = -9.457425049, c = 2.635173636,  d = -3.098230243}

Is Maple able to find a formal solution of this same sys(1) system?

solve(sys(1), unknowns);

I waited about 5 minutes (Maple 2015) and couldn't get one.
So my guess is that Maple can't find a solution in terms of (here) x (alone)
 

Can you give the expression of the general term of this serie?
That is, if one write u(x, t)=u0(x) + u1(x)*t + ...un(x)*t^n+ ..., what is the expression of un(x)?
Or if you prefer, what is the recurrence relation between un(x) and un-1(x), un-2(x), ...un-p(x), (n>p)?

If we take the problem from the other side, one observe that the taylor expansion of u(x,t) with respect to t has its first terms equal to the first ones of the serie you give.

u := 1/(1+exp(x-5*t))^2:
S := map(simplify, taylor(u, t))

Of course it's easy do to what I'm going to do once we accept the previous result.

  • the coefficients of the successive terms are 10, 50/2, 250/6, 1250/24, ...
    thus terms of the form 10 * 5^k / k! with k=0..+oo
     
  • this suggest that the serie could be a Taylor expansion of some function w(x, t)
    moreover, giving the expression of the terms that contain t, this Taylor expansion, which should contain terms like (t-a)^k can be done at a=0
     
  • for k=0 one has w(x, 0) = f(x) = 1/(1+exp(x))^2.
    let us observe a few derivatives of f(x) with respect to x:
    f := (1/(1+exp(x))^2:
    print~([seq(simplify(eval(diff(f, [x$k]), t=0)), k=0..5)]):
    
    we observe that these derivative are equal, up to a miltiplicative constant, to the terms that contain x in your serie
     
  • if w(x, t) is of the form w(x+c*t), it's easy to see that 
    # set z=x+c*t
    diff(w(x+c*t), t$k) = diff(w(z), z$k)*c^k
    this could lead to the inference that your serie is the Taylor expansion (1/(1+exp(x+c*t))^2 around t=0
     
  • under this hypothesis one easily finds that c=5


If the assumption that your serie is such a Taylor expansion, it converges for any real values of x and t.

@Michael_Watson 

restart:

e := L/(z*sqrt(z^2 + L^2))

 

# case L >> z
# set L=a*z with a >> 1

eval(e, L=a*z);
numer(%)/z / simplify(denom(%)/z) assuming z > 0;
eval(%, a=L/z);

asympt(%, L, 2):
expand(simplify(%)) assuming z > 0

 

L/(z^2*(L^2/z^2+1)^(1/2))

 

1/z+O(1/L^2)

(2)

# case z >> L
# set z=a*L with a >> 1

eval(e, z=a*L);
numer(%)/L / simplify(denom(%)/L) assuming L > 0;
eval(%, a=z/L);
asympt(%, z, 4):
expand(simplify(%)) assuming L > 0

 

1/(z*(z^2/L^2+1)^(1/2))

 

L/z^2+O(1/z^4)

(3)

 


 

Download asympt.mw

One way with geom3d package

EDITED VERSION

 

restart:
with(geom3d):

# first line

plane(P1, 16*x-2*y-11*z, [x,y,z]):
plane(P2, 14*x-y-10*z-3, [x,y,z]):
line(L1, [P1,P2]):
Param_Eq__1 := Equation(L1,'t');
                    [1                     ]
                    [- + 9 t, 4 + 6 t, 12 t]
                    [2                     ]
# second line (=L2)
# let u be the common value of (x-2)/3, (y-5)/2 and (z-2)/4

X := solve((x-2)/3=u, x);
Y := solve((y-5)/2=u, y);
Z := solve((z-2)/4=u, z);
r := u=solve(X=rhs(Param_Eq__1[1]), u);
eval([X, Y, Z], r)
                               1      
                         u = - - + 3 t
                               2      
              [1                     ]
              [- + 9 t, 4 + 6 t, 12 t]
              [2                     ]

Thus lines L1 and L2 are  identical.

First problem only (I'm going to bed)

Step 1: define M and N this way

restart:
M := (a, b) -> (a+b)/2;
N := (a, b) -> a*b/M(a, b);

 

Step 2: the recurrence is

a(n+1) = M(a(n), b(n)),  b(n+1) = N(a(n), b(n))

Set  z(n+1) = a(n+1) * b(n+1) and observe that z is invariant under the transformations M and N.

z(n+1) = a(n+1) *  b(n+1) = a(n) * b(n) = ... = a(0) * b(0)



Step 3 : find a limit (aka stationary) solution by solving this system of equations (could be done by hand)

solve({M(x,y)=x, N(x, y)=y})
                         {x = y, y = y}

Thus a stationary solution is such that x=y (or limit(a(n), n=+oo) =  limit(b(n), n=+oo) if you prefer


Step 4 : inject this stationary solution within z : 
As limit(a(n), n=+oo) =  limit(b(n), n=+oo), one gets limit(z(n), n=+oo) = z(0) = a(0)*b(0) and thus 

limit(a(n), n=+infinity) = limit(b(n), n=+infinity) = sqrt(a(0)*b(0))


Example : a(0)=1 and b(0)=2 gives  limit(a(n), n=+oo) =  limit(b(n), n=+oo) = sqrt(2)

 

 

Does this suit you?


(u instead of u(y))

restart
alias(u=u(y, z));
de := diff(u, y$2)+3*delta*diff(u, y)^2*diff(u, y$2)=p;

# and for subsequent computations
eval(de, u=U(y));
# example
dsolve(%, U(y));

Note that the trick alias(u(y)=u(y, z)) won't produce the desired result

First 23 24 25 26 27 28 29 Last Page 25 of 47