## 17 Badges

7 years, 115 days

## @MaPal93Neither your "simple" ...

@MaPal93

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
230323_updated1_Problem_mmcdara_simple.mw

## @MaPal93 Let me be clear: You asked...

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.

## @MaPal93 What do you want to achiev...

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

`sqrt(Sigma__0jk)*(sqrt(Sigma__0ki)*beta__jk3*rho__XX__23+sqrt(Sigma__0jk)*beta__jk2)/(2*sqrt(Sigma__0jk)*sqrt(Sigma__0ki)*beta__jk2*beta__jk3*rho__XX__23+Sigma__0jk*beta__jk2^2+Sigma__0ki*beta__jk3^2+sigma__ujk^2)`

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.

UPDATED VERSION
230323_updated1_Problem_mmcdara.mw

## @rcorless I suppose that a post abo...

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.

## @MaPal93 I will look to your file t...

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.

 > restart
 > with(Statistics):
 > # 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);           b*thisproc(a)      else 'procname'(e)      end if end proc: #------------------------------------------------------------------------ TypeTools:-AddType(      RV,      {RandomVariable,      'RandomVariable^posint',      '`*`'({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:     else       a := eval(a, E = (RV -> Mean(RV))):     end if;   else     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       else         a := eval(t, E = (RV -> Mean(RV))):       end if;       b := 1     else       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);
 (1)
 > u := RandomVariable(Abstract(0, sigma__u^2)); (Mean, Variance)(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);
 (3)
 > # Unfortunately Maple can compute this Mean(R);
 (4)
 > # 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) print(): FOC := mu = rhs(%);
 (5)
 > # 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));
 (6)
 > #   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)]);
 (7)
 > # 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)
 (8)
 > # Find the value of lambda which minimizes < VAR > # # Relation (8b) in the paper solve(diff(VAR, lambda), lambda): Lrel := 'lambda' = %;
 (9)
 > # Relation (9a) in the paper eval(FOC, [alpha=-mu/2/lambda, beta=1/2/lambda]): murel := isolate(%, mu)
 (10)
 > # How do the authors get their relation (9b) ???? eval(FOC, murel): isolate(%, lambda); eval(Lrel, beta=1/2/lambda): isolate(%, lambda);
 (11)
 >

Download Abstract_RV.mw

## @MaPal93 "What this equation s...

"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.

## What do you want to achieve in a few wo...

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.

 > restart
 > with(Statistics):
 > 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;
 (1)
 > first_order_condition := mu = collect(Mean(R), p__0);
 (2)
 > var := Variance(R)
 (3)
 > var := eval(var, sigma__0^2= Sigma__0)
 (4)
 > L := solve(diff(var, lambda)=0, lambda): 'lambda' = L;
 (5)
 > isolate(eval(first_order_condition, alpha=-mu/2/lambda), mu); isolate(eval(%, beta=1/lambda), mu)
 (6)
 >

Download What_I_unerstand.mw

## @rcorless Hi, this could be ve...

@rcorless

Hi,
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...

@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

 > restart:

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(                      eval(                        convert(                          taylor(op(0, f)(op(1, f)), op(1, f)=A, m),                          Diff                        ),                        op(1, f)=A+s*H),                      polynom                    ):   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)])[];   else     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"):     return   else     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:
 > Describe(CDDF)
 # 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);
 (1)
 > convert(simplify( convert(taylor(c, h=0), polynom) ), Diff);
 (2)
 >

Download NotSimplified.mw

## @jetboo  About the double arro...

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

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

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*~p;
3*p;
[3, 6]
[3, 6]
s*p;
s*~p;
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).

formulas_example.mw

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

Download CenteredDividedDifferences.mw

## @dharr Thanks for your comments.(An...

@dharr

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.

## Building all the cycles...

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

 > restart:
 > with(GraphTheory):
 > CG := SpecialGraphs:-CageGraph(3,5): DrawGraph(CG);
 > cb     := CycleBasis(CG); cycles := map(x->Edges(InducedSubgraph(CG, x)), cb):
 (1)
 > symmdiff(cycles[[1, 2, 3, 5]][]);
 (2)
 > # the result contains two disjoint cycles [1, 2, 9, 8, 5] , [3, 4, 10, 6, 7]

Download CounterExample.mw

## @Rouben Rostamian   Let's wait...

@Rouben Rostamian

Let's wait for some clarification.

## @Rouben Rostamian  @Rouben Rostamia...

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 I've never seen a pr...

@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).

 > restart:
 > with(GraphTheory):
 > 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)
 > # 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);
 (2)
 > 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)[]]
 (3)
 > # Check GG := Graph({convert~(A1, set)[]}): DrawGraph(GG);
 >

Download CycleBasis_2.mw

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
﻿