5269 Reputation

17 Badges

7 years, 115 days

MaplePrimes Activity

These are replies submitted by mmcdara


Neither your "simple" problem nor your "complex" problem has a solution.

Given the content of your last message, I understand that you havn't even read the last file I sent you.
Had you, you would have understood that the "simple" version does not have a solution.

So here is this same file with a name that will suit you more


Let me be clear: You asked a question and I gave a reply, then you modified your initial question and gave another reply, and you keep modifying the question and keep asking for help.
I have already spent several hours answering your successive answers without you giving me any retribution like "I vote up", or "I select this answer as the best one".

So I won't try to solve this umpteenth variant, given that your previous problem doesn't have any solution.


What do you want to achieve with this command:

eqmujk := eval(mu__jk, %[]); # ?

It is not syntactically correct for %[ ] means "the elements of the previous output", which assumes this previous output i euther a list or a set, which is not the case.
Here is the output %[ ] refers to


I suppose you have just reproduce the syntax which produces output on line (17), but here  %[ ] refered to sol1 which was a list and "extracted" the inner list 

[mu__jk = ..., lambda__jk = ...]

Then the command 

eqljk := eval(lambda__jk, %[]);

meant "evaluate lambda__jk considering the the rules [mu__jk = ..., lambda__jk = ...].

I could have written instead, and I realize now that my initial syntax was confusing for newbies (no effence intended).

eqljk := rhs(sol1[][2])

If I had written the things that way you would probably defined  eqmujk as

eqmujk := rhs(sol1[][1])

I am right?

Further down in the text : 

eqlki := rhs(sol1[][2]);
eqmuki := rhs(sol2[][2]);

This fixes several errors.
I send you a corrected file, but the last solve command, which I rewrote correctly, doesn'tt provide any solution.



I suppose that a post about the subject will be highly appreciated.

I always used modified equations for PDE. 
Thirty years ago my area of interest was fluid mechanics, more precisely the design of numerical schemes for Large Eddy Simulation, capable to capture the Kolmogorov's energy cascade in Homogenous Isotropic Turbulence.
The modified equation was the main tool to discard classical toolsin the domain, to promote some others and to imagine new one.
But a lot of water has passed under the bridge since then.


I will look to your file tomorrow and keep you informed if I have some ideas or anu difficulty to understand it.

Meanwhile here is a more detailed version of my previous worksheet.



# Definition of < E > (Expectation) and of type < RV > has been
# provided by Carl Love

E:= proc(e::algebraic)
     local a,b;
     if not hastype(e, RV) then e
     elif e::RV then 'procname'(e)
     elif e::`+` then map(thisproc, e)
     elif e::`*` then
          (a,b):= selectremove(hastype, e, RV);
     else 'procname'(e)
     end if
end proc:


     '`*`'({RandomVariable, 'RandomVariable^posint'})

# This procedure is aimed to assess the variance of a linear combination
# of random variables of < Abstract > distitribution

reduce := proc(t)
  local a, b:
  if type(t, `*`) then
    a, b := selectremove(has, t, E):
    if type(op(1, a), `^`) then
      a := eval(a, E=(RV -> (Variance+Mean^2)(op(1, RV))))
    elif type(op(1, a), `*`) then
      a := 0:
      a := eval(a, E = (RV -> Mean(RV))):
    end if;
    if has(t, E) then
      if type(op(1, t), `^`) then
        a := eval(t, E=(RV -> (Variance+Mean^2)(op(1, RV))))
      elif type(op(1, a), `*`) then
        a := 0
        a := eval(t, E = (RV -> Mean(RV))):
      end if;
      b := 1
      a := 1:
      b := t:
    end if:
  end if:
  return a*b
end proc:

Abstract := (m, v) -> Distribution(Mean=m, Variance=v);
v := RandomVariable(Abstract(p__0, Sigma__0));

(Mean, Variance)(v);

proc (m, v) options operator, arrow; Statistics:-Distribution(Statistics:-Mean = m, Statistics:-Variance = v) end proc




p__0, Sigma__0


u := RandomVariable(Abstract(0, sigma__u^2));

(Mean, Variance)(u);



0, sigma__u^2


# Define a random variable R as a linear combination of u and v.
# Relation (8a) in the paper

R := (1-lambda*beta)*v - lambda*(u +alpha);



# Unfortunately Maple can compute this




# So do the things this way:
#   First : apply the linear oprrator < E > to < R >.
#   Next  : replace < E > by < Mean > in the previous result.

eval(E(R), E = (RV -> Mean(RV))):
'E(R)' = collect(%, p__0);

# Here is what you call FOC (First Order Condition)

FOC := mu = rhs(%);

E(R) = (-beta*lambda+1)*p__0-lambda*alpha



mu = (-beta*lambda+1)*p__0-lambda*alpha


# Proceed the same way as befor to compute the second non-centered
# moment < E(R^2) > of R
#   First : apply the linear oprrator < E > to < R^2 >.

E2 := E(expand(R^2));



#   Next  : the evaluation of this expression is a little bit tricky and
#           requires using the procedure < reduce >equation (9b) ????

E2 := add(reduce(x), x in [op(E2)]);



# To get the variance of R, just write that
#          Var(R) = E((R-E(R))^2) = E(R^2) - (E(R))^2

VAR := simplify(E2 - rhs(FOC)^2)



# Find the value of lambda which minimizes < VAR >
# Relation (8b) in the paper

solve(diff(VAR, lambda), lambda):
Lrel := 'lambda' = %;

lambda = beta*Sigma__0/(beta^2*Sigma__0+sigma__u^2)


# Relation (9a) in the paper

eval(FOC, [alpha=-mu/2/lambda, beta=1/2/lambda]):
murel := isolate(%, mu)

mu = p__0


# How do the authors get their relation (9b) ????

eval(FOC, murel):
isolate(%, lambda);

eval(Lrel, beta=1/2/lambda):
isolate(%, lambda);

p__0 = (-beta*lambda+1)*p__0-lambda*alpha


lambda = 0


lambda = 0





"What this equation seems to represent." Why are you introducing a variance here? I want to minimize the expectation of something in the form (a-b-c)^2, where (a-b-c) includes RVs and constants. 

Line 1 of the paper is 

E( (X-mu)^2);

# with 

X = (1-lambda*beta)*v - lambda*(u+alpha) 

# and it happens, equation (8a) that 

mu = (1-lambda*beta)*p0 -lambda*alpha

# which is exactly

mu =  (1-lambda*beta) * E(v) - lambda * E(u+alpha)

# thus mu=E(X) and 

E( (X-mu)^2) = E( (X-E(X))^2) - = (bydefinition) Variance(X)
  • "Equation (7) is just the expansion of equation (6) given the expectationsn variances and correlation oof v tilde and u tilde."

Yes, but given that the correlation between vtilde and utilde is 0 (by definition)

  • I understand that what is called "first order condition" is nothing but the expectation ("mean") of"

No. First order conditions (FOCs) here mean partial derivatives set to 0. For example, FOC w

You're right about the general definition of FOC...
... but this doesn't change anithing to the fact that E( (X-c)^2) (X being a random variable and c a number) is minimal when c = E(X).
This is a mathematical theorem we encounter in several different forms in other fields than Probability.
For instance, in Mechanics, the inertia of a solid around a point is minimal vhen this point is the center of mass of the solid.
It seems to me that the paper you reproduce looks for complications.

In the excerpt of the paper you provide:

  • E seems to denote the "Mathematical Expectation" and E[Q] seems to represent the expectation of the random variable Q?
  • Here Q is a combination of 2 random variables v tilde and u tilde , whose expectations and variances are given on the line which follows equation (6).
    What this equation seems to represent
    Variance(v*(1-lambda*beta) - lambda*u - lambda*alpha)
  • Equation (7) is just the expansion of equation (6) given the expectationsn variances and correlation oof v tilde and u tilde.
  • So youj search for the value of lambda such that the variance of RV is minimal.
  • I understand that what is called "first order condition" is nothing but the expectation ("mean") of
    v*(1-lambda*beta) - lambda*u - lambda*alpha

So please, if you want some help explain claitly what you want.
If it is working with random variables it's likely one can do the things in a simpler way.

Here is a code which details what I said above.
I'm quite sure it is exactly what you want, but it might help you formalize the problem in a better wat.



v := RandomVariable(Normal(p__0, sigma__0)):
u := RandomVariable(Normal(0, sigma__u)):

# Is the next line right?

R := v*(1-lambda*beta)-lambda*u - lambda*alpha;



first_order_condition := mu = collect(Mean(R), p__0);

mu = (-beta*lambda+1)*p__0-lambda*alpha


var := Variance(R)



var := eval(var, sigma__0^2= Sigma__0)



L := solve(diff(var, lambda)=0, lambda):
'lambda' = L;

lambda = beta*Sigma__0/(beta^2*Sigma__0+sigma__u^2)


isolate(eval(first_order_condition, alpha=-mu/2/lambda), mu);
isolate(eval(%, beta=1/lambda), mu)

mu = 2*(-beta*lambda+1)*p__0


mu = 0





this could be very interesting.

About 25 years ago I was working on a specific topic named "equivalent equation" (sometimes named "modified equation").
I don't know if you are familiar with this?

In a few words, let A a differential operator (linear to simplify), u some function, and A(u)=0 a partial differential equation.
Let a numerical scheme, an approximation of  u, and S(v)=0 an approcumation of A(u)=0 .
The truncation error is a classical characterization of S, this is basically what the mtaylor function in my post provides after having substracted A(u).

The "equivalent equation" searches what continuous operator Aeq the numerical scheme S approximates with a null truncation error.
This is a kind of inverse problem whose solution is very helpful to understand how a numerical scheme S behaves the way it does (for instance to prove it diffusive or anti-diffusive).

I did some work with Maple ( Maple V I think :-) ) and went back to it 10 years ago.
I still have some material about this if you are interested.

Thanks for your comment.

@Preben Alsholm 

I had exactly the same problem with Maple 2015.2 while using taylor, which is why I use mtaylor instead.

I identified that taylor returning an expression with O(...) was the cause (don't ask my why) for Maple not simplifying the formulas.
Converting into polynom before converting onto Diff/diff fix this issue, but I nevertheless kept mtaylor in my post.

Look here


CDDF stands for Cendered Divided Difference Formula

CDDF := proc(f, A, H, n, stencil)
  description "f = target function,\nA = point where the derivatives are approximated,\nH = step,\nn = order of the derivative,\nstencil = list of points for the divided differenceCDDF\n";
  local tay, p, T, sol, unknown, Unknown, expr:

  tay := (s, m) -> convert(
                         taylor(op(0, f)(op(1, f)), op(1, f)=A, m),
                       op(1, f)=A+s*H),

  p   := numelems(stencil):
  T   := add(alpha[i]*tay(i, p+1), i in stencil):
  T   := convert(%, diff):

  if p > n+1 then
    sol := solve([seq(coeff(T, h, i)=0, i in subsop(n+1=NULL, [$0..p]))], [seq(alpha[i], i in stencil)])[];
    sol := solve([seq(coeff(T, H, i)=0, i in subsop(n+1=NULL, [$0..n]))], [seq(alpha[i], i in stencil)])[];
  end if:

  if `and`(is~(rhs~(sol)=~0)[]) then
    WARNING("no solution found"):
    unknown := `union`(indets~(rhs~(sol))[])[];
    Unknown := simplify(solve(eval(T, sol) = Diff(op(0, f)(A), A$n), unknown));
    sol     := lhs~(sol) =~ eval(rhs~(sol), unknown=Unknown);
    expr    := normal(eval(add(alpha[i]*op(0, f)(A+i*H), i in stencil), sol));
  end if:

  return expr
end proc:


# f = target function,
# A = point where the derivatives are approximated,
# H =
# step,
# n = order of the derivative,
# stencil = list of points for the divided
# differenceCDDF
CDDF( f, A, H, n, stencil )


# 2-point approximation of diff(f(x), x) | x=a

c := CDDF(f(x), a, h, 1, [-1, 1]):

simplify(taylor(c, h=0));
convert(simplify(taylor(c, h=0)), Diff);
simplify(convert(taylor(c, h=0), Diff));
convert(simplify(mtaylor(c, h=0, 4)), Diff);



series(Diff(f(a), a)+((1/6)*(Diff(Diff(Diff(f(a), a), a), a)))*h^2+((1/120)*(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a)))*h^4+O(h^5),h,5)


series(Diff(f(a), a)+((1/6)*(Diff(Diff(Diff(f(a), a), a), a)))*h^2+((1/120)*(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a)))*h^4+O(h^5),h,5)


Diff(f(a), a)+(1/6)*(Diff(Diff(Diff(f(a), a), a), a))*h^2


convert(simplify( convert(taylor(c, h=0), polynom) ), Diff);

Diff(f(a), a)+(1/6)*(Diff(Diff(Diff(f(a), a), a), a))*h^2+(1/120)*(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a))*h^4







About the double arrow procedure: f and g above are two different ways to do the same thing

f := (x, n) -> x^n:
f(u, 3)
g := n -> x -> x^n:

So, why did I choose to  to write 

tau := (d, s) -> p -> p + s*~[2-d, d-1]

instead of

tau := (d, s, p) -> p + s*~[2-d, d-1]


Because, and this is questionable, I felt that the double arrow definition matched better the idea to "apply a translation along direction d of amount s to the point p" than a single arrow definition would do.

About the < *~ > operation: this is a "component wise operation", and it means "apply the multiplication component by component".
If s was numeric the  < *~ > would be useless and could be replaced by  a simple < * >. But this is not the case if s is not numeric:

p := [1, 2]:

                             [3, 6]
                             [3, 6]
                            s [1, 2]
                            [s, 2 s]

In case you would be interested, the attached file, which I tried to be pedagogical, explains how to construct divided difference formulas and to assess their order of approximation:

  • The first example builds a 2-point forward divided difference formula to approximate diff(f(x), x).
  • The first example builds a 3-point centereddivided difference formula to approximate diff(f(x), x$2).
  • The third example builds a 4-point centered divided difference formula to approximate diff(f(x), x$3).
  • The third example builds a 5-point centered divided difference formula to approximate diff(f(x), x$3).

The file below contains a procedure which builds cenered divided difference formulas in 1D



Thanks for your comments.

BTW, I do remember the work you did on the blog you mention, but I hadn't tied it to this discussion.
I'll take another look at it.

@Carl Love 


At @dharr : your "combining cycles" solution is very astute and works well on the case I proposed, but here is a counter example



CG := SpecialGraphs:-CageGraph(3,5):


cb     := CycleBasis(CG);
cycles := map(x->Edges(InducedSubgraph(CG, x)), cb):

[[1, 2, 3, 4, 5], [1, 2, 3, 7, 6], [1, 5, 4, 10, 6], [1, 5, 8, 7, 6], [1, 2, 9, 8, 5], [1, 2, 9, 10, 6]]


symmdiff(cycles[[1, 2, 3, 5]][]);

{{1, 2}, {1, 5}, {2, 9}, {3, 4}, {3, 7}, {4, 10}, {5, 8}, {6, 7}, {6, 10}, {8, 9}}


# the result contains two disjoint cycles

[1, 2, 9, 8, 5] , [3, 4, 10, 6, 7]



@Rouben Rostamian  

Let's wait for some clarification.

@Rouben Rostamian
@Madhukesh J K

You right Rouben. 

I began changing CreateVector into Vector and CreateMatrix into Matrix, but the OP uses indices in the range 0..N.
So I changed that into 

K := Array(0 .. N, 0 .. N, 0);
F := Array(0 .. N, 0);

But there are things I still don't understand:

  • Where is the differential equation?
  • The integral formula (8) can't be a weak formulation for v(x), which I suppose is the projection function appears only in one integral.
  • Is it a Galerkin method with the projection function of the same form of those in the ansatz of u(x)?
  • As the problem involves (I suppose a fourth order derivative of u(x)... maybe a non linear problem of Mechanics?) the test functions can't be linear.
  • What does the procedure Llinear do? 

@Carl Love 

Thanks for you comments.

I've never seen a proper formal definition of cycle basis  neither did I.
Of course all the cycles in a graph can be constructed through elementary operations (which doen't mean it's easy to code this).

Do you know how to get the cycle [3, 4, 5, 6, 8] for your basis?
This is the kind of problem I often encountered and was never able to solve in a generic and robust way (I always found a test case for which my code did not give the right answer).

For the case described in my previous file here is the (quite ad hoc) way to build  the cycle [3, 4, 5, 6, 8] (unfortunately, I'm off work for two weeks and I can't get you the codes I wrote to build other cycles from the CycleBasis).



C1 := {seq({i, i+1}, i=1..6), {7, 1}}:
C2 := {{1, 8}, {8, 3}, {8, 6}}:
G := Graph(C1 union C2):
DrawGraph(G, style=planar);

CB := CycleBasis(G);


[[1, 2, 3, 8], [1, 2, 3, 4, 5, 6, 7], [1, 7, 6, 8]]


# Start from CB[2] and try to obtain the cycle [5, 6, 8, 3, 4]
# CB[2] contains the subsequence [1, 2, 3] which is also in CB[1].
# So we can replace 1, 2, 3 in CB[2] by the other path from 1 to 3 in C1
# that is 1, 2, 3 by 1, 8, 3.
# To do this, delete the edges {1, 2} and {2, 3}

G1 := Graph(Edges(G) minus { {1, 2}, {2, 3}});
DrawGraph(G1, style=planar);

CB1 := CycleBasis(G1);


GRAPHLN(undirected, unweighted, [1, 3, 4, 5, 6, 7, 8], Array(%id = 18446744078204485750), `GRAPHLN/table/64`, 0)



[[1, 7, 6, 5, 4, 3, 8], [1, 7, 6, 8]]


L1 := ListTools:-Rotate(CB1[1], -1);
L2 := ListTools:-Rotate(CB1[2], -1);

# L2 is a sublist of L1
# Step 1: build the arcs corresponding to L1:

A1 := [ seq([L1[i], L1[i+1]], i=1..nops(L1)-1), [L1[-1], L1[1]] ];

# Step 2: build the arcs corresponding to L2:

A2 := [ seq([L2[i], L2[i+1]], i=1..nops(L2)-1), [L2[-1], L2[1]] ];

# Step 3: remove in A1 all the arcs which are in A2

R2 := copy(A2):
for a in A2 do
  p := ListTools:-Search(a, A1):
  if p <> 0 then
    A1 := subsop(p=NULL, A1):
    R2 := subsop(1=NULL, R2):
  end if
end do;

A1, R2;

# Step 4: add to A1 the reverse arcs that are still in R2

A1 := [A1[], ListTools:-Reverse~(R2)[]]

[8, 1, 7, 6, 5, 4, 3]


[8, 1, 7, 6]


[[8, 1], [1, 7], [7, 6], [6, 5], [5, 4], [4, 3], [3, 8]]


[[8, 1], [1, 7], [7, 6], [6, 8]]










[[6, 5], [5, 4], [4, 3], [3, 8]], [[6, 8]]


[[6, 5], [5, 4], [4, 3], [3, 8], [8, 6]]


# Check

GG := Graph({convert~(A1, set)[]}):





REMARK: A simpler way would be to:

  1. remove the arcs [6, 7] and [7, 1] the same way I used to build graph G1,
  2. apply the same procedure to remove arcs [1, 2] and [2, 3].

In this case, we use the same procedure twice, but this can only be done because we have seen what the structure of the graph was: this ia a tailored procedure, not a generic one. 

First 12 13 14 15 16 17 18 Last Page 14 of 113