Carl Love

## 24663 Reputation

10 years, 56 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

## Initial condition...

I can't say for sure without seeing it in context, but I suspect that N1 is an initial condition, for example, the value of the formula for n=1. Every recursive formula will have at least one initial condition, although the value of n where it occurs isn't necessarily 1.

The Maple command rsolve can solve some recursion formulas that can be expressed via algebraic equations (in which case they're usually called recurrence relations or difference equatios). Observe here how the specification of initial conditions affects the results:

```#No initial condition specified, so Maple chooses N(0) as the initial condition:
rsolve({N(n) = a+N(n-1)}, N(n));
N(0) - a + a (n + 1)

#Initial condition specified at n=0:
rsolve({N(n) = a+N(n-1), N(0)=N0}, N(n));
N0 - a + a (n + 1)

#Initial condition specified at n=1:
rsolve({N(n) = a+N(n-1), N(1)=N1}, N(n));
N1 + a (n + 1) - 2 a

```

## LinearSolve(A, B) instead of A^(-1).B...

I suspect that the memory explosion occurs when you do the matrix inverse (KL + a0.ML + a1.CL)^(-1). A matrix multiplication of the form X:= A^(-1).B is mathematically equivalent to X:= LinearAlgebra:-LinearSolve(A, B), but the latter is usually much more efficient.

## Tweaked...

I think that this does what you asked; let me know:

```IsLabeledIsomorphicSubgraph:= proc(G1::Graph, G2::Graph)
uses GT= GraphTheory;
local V:= GT:-Vertices(G1), R;
if {V[]} subset {GT:-Vertices(G2)[]} then
R:= GT:-IsSubgraphIsomorphic(G1, GT:-InducedSubgraph(G2, V), 'isomorphism');
if [R][1] then subs(R[2], GT:-Edges(G1)) else false fi
fi;
false
end proc
:```

## Almost certainly inconsistent...

• Are you implying that when Maple's solve() seems to get stuck is because the system has no solution?

No, I'm implying that when there are more equations than unknowns, even just one more, then the system is almost certainly inconsistent, unless you have some strong theoretical reason for believing otherwise. Imagine drawing 3 curves in a plane. It's almost certain that there won't be a single point where all three intersect.

If we can find any inconsistent subset, under any values of parameters, then the whole system is inconsistent. So, the plan below is

1. Let Eq be the set of 12 equations (assumed here to be expressions implicitly equated to 0); let be the set of 6 variables that you want to solve for.
2. Extract the numerators of the equations; the denominators are irrelevant. Call this EqN.
3. Set any variables in EqN that are not in to 1, or any other simple, exact value, possibly 0. Call this EqP.
4. Find subsets of EqP of size 6 that contain all variables in V. Call this Eq6.

The above steps are very easy. They'll take 5 seconds or less in total, like this:

EqN:= numer~(Eq):
EqP:= eval(Eq, indets(EqN, name)=~ 1):
Eq6:= select(E-> V subset indets(E, name), combinat:-choose(EqP, 6)):

5. Take any member of Eq6 and try solve on it:

SolE:= {solve}(Eq6[1], V);

6. Remove any solution that contains free variables:

SolE2:= remove(s-> ormap(evalb, s), SolE);

7. If SolE2 is empty, go back to step 5 and pick another set of 6.
8. Take any solution in SolE2 and use it to evaluate the other 6 equations:

EqV:= simplify(eval(EqP minus Eq6[1], SolE2[1]));

9. If EqV = {0}, then, miraculously, the system is consistent, and you have a solution, at least for the parameter valuation used in step 3.
10. If EqV contains anything that is obviously not equal to 0, then the entire system is inconsistent; you're done; period.

• Did you check my script?

Yes.

## unapply for indirect parameter reference...

I call something like f:= x*cos(x); ff:= x-> f, an attempted indirect parameter reference, the x being the parameter. Those references don't work for procedures created with -> or proc. The most-common way to handle them is ff:= unapply(f, [x]). To be more robust, you should pass x into q.

The quotes that you put around x don't do anything.

## expand, coeffs, log, simplify, sort...

It can be done like this:

#Important! The next line has no minus sign!
C:= coeffs(expand(lhs(pde)), e^(gamma*tau), 'n'):
sort(`[]`~(simplify(log[e^(-gamma*tau)]~([n]), symbolic), [C]));

The returned expression is a list of pairs. The 1st element of each pair is the power to which e^(-gamma*tau) is raised; the 2nd is the corresponding coefficient.

## Without strings or parse...

When parse returns an expression containing variables (i.e., names, identifiers, function names), and those variables weren't originally global, they won't be the same as the ones passed in.

The two types of substitution that you want can also be done with a simple loop, of course:

```SublistSubs:= proc(S::(list=anything), L::list, {type2::truefalse:= false})
local s:= lhs(S), t:= rhs(S), r:= nops(s), T:= Vector(r--, [t], fill= ()), k, V:= <L>;
for k to nops(L)-r do
if L[k..k+r] = s then
if type2 then V[k]:= t; k+= r else V[k..(k+= r)]:= T fi
fi
od;
[seq](V)
end proc
:
map[3](
SublistSubs,
[a,b,c]= x, [2,4,5,1,a,4,5,a,b,434,a,b,c,5,5,5,1],
type2=~ [true,false]
);
[[2, 4, 5, 1, a, 4, 5, a, b, 434, x, 5, 5, 5, 1],
[2, 4, 5, 1, a, 4, 5, a, b, 434, x, b, c, 5, 5, 5, 1]]

```

## Simplify in exact arithmetic...

For any values of kB, and R, both x=-1 and x=1 are roots. No decimal approximations are needed to verify this; it can be done in exact arithmetic with symbolic k, B, and R:

simplify(eval~(P, x=~ [-1,1]));
[0, 0]

Your expression initially appears quite lengthy. That's superficial. It can be vastly reduced in size by

P1:= simplify(normal(P));

And for the purpose of understanding its roots, it can be further reduced by

P1n:= numer(P1);

The expression is now easy to understand, and you can see why those are the roots for x. It's not totally obvious, but could easily be verified by, for example, pencil and paper and only high-school algebra.

## Prime testing is much faster than factor...

It is surprising---but nonetheless true---that determining whether or not a large number is prime can be done in much less time than factoring a number of similar size that has at least two large prime factors. That fact is fundamental to the usefulness of the RSA encryption method. So the fact in the first sentence is fairly well known because that encryption method is of immense practical importance.

Let's see just how fast it is. First, let's estimate the number of 116-digit primes:

evalf(Int(1/ln(x), x= 1e115..1e116));
3.37895657556787*10^113

That's far more than the number of atoms in the Universe or even the volume of the observable Universe in cubic nanometers. So, if we were to choose one of these numbers at random, the chance that it could've been pre-stored is infinitesimal, right?  Surely it'd take more than one atom to store a 116-digit number.

```R:= rand(10^115..10^116-1):
randprime:= proc()
local n, p:= R();
for n while not isprime(p) do p:= R() od;
n, p
end proc
:
(n,p):= CodeTools:-Usage(randprime());
memory used=0.71MiB, alloc change=0 bytes,
cpu time=0ns, real time=5.00ms, gc time=0ns

n, p := 210, 24729845748629285111963221649430106890373391188328886700762826051452963665879923038495438270815820598184513197776291
```

So, it took 5 milliseconds to generate and test 210 random 116-digit numbers, the last of which was the first prime among them. So generating and testing 116-digit numbers takes about 25 microseconds apiece.

## Need to solve the 2 ODEs together...

It's possible to solve your ode1 by itself because it doesn't contain theta(x). But it's not possible to solve ode2 by itself because it contains 2 unknown functions, f(x) and theta(x). Nor is there any need to solve it by itself; just solve them together:

b5:= subs(a5, {ode1, ode2});
c5:= dsolve({b5[], bcs, bcs1}, numeric);

And do likewise for b6 and c6.

## Using matrices...

By using matrix multiplication and matrix stacking, the code can be shortened to this:

```restart:
Digits := 15:
T5Scheme:= proc(XY0::Matrix, a::Matrix, L::nonnegint)
local i, k, N:= upperbound(XY0)[1], XY:= rtable(XY0);
for k to L do XY:= <seq(a.XY[i-2..i+2], i= 3..3^(k-1)*(N-6)+4)> od;
XY
end proc
:
## Initial polygon
q2:= -sqrt(2): q3:= sqrt(3):
XY:= <
-2,  1, 2,  0, q2, -2,  1, 2,  0, q2, -2,  1, 2,  0, q2  | #x
0, q3, 0, -2, q2,  0, q3, 0, -2, q2,  0, q3, 0, -2, q2    #y
>/2:
st2:= <-2, 1, 2, 0, q2, -2  |  0, q3, 0, -2, q2, 0>/2:
L:= [13, -22, -81, -107, 179, 729, 1233]/1296:
a:= Array(-6..8, [L[], ListTools:-Reverse(L)[], 0]):
J:= [seq](6..-6, -3):
A:= <a[J], a[J+~1], a[J+~2]>:
st1:= T5Scheme(XY, A, 3):
plot([st1, st2], linestyle= [1, 3], color= [black, red]);
```

## Use the numerators...

In almost all Maple commands, using an expression that's not an equation as an equation is equivalent to setting that expression equal to 0. Solving an expression set equal to 0 is almost equivalent to just solving its numerator set to 0, the only difference being that you might want to validate the solutions by checking whether they make the denominator 0, possibly under some restrictions. So, this works for your case:

pdsolve(numer~([eq1, eq2]));

Regarding your literal question "What does this error mean?": I don't know. It's very often much easier for me to correct or prevent an error in Maple than to understand it. It's a good question, but answering it might take several hours of investigating, digging through the code.

## Procedure for any of the 2^8 graph produ...

Using the same 1--8 numbering scheme as @lcz posted (from  Barik, S., Bapat, R.B., Pati, S.: "On the Laplacian spectra of product graphs" Appl. Anal. Discrete Math. 9, 39–58 (2015)), here is a procedure for any of the 2^8 = 256 possible products of two graphs:

```RelationGraph:= proc(f, V::list({integer, symbol, string}))
local n:= {\$1..nops(V)}, i, F:= i-> j-> f(V[i],V[j]);
GraphTheory:-Graph(V, [seq](select(F(i), n, _rest), i= n))
end proc
:
GraphProducts:= proc(
G::Graph, H::Graph,
Rel::{
set(integer[1..8]),
identical(
"Cartesian", "box", "strong", "normal", "AND", "modular", "homomorphic",
"lexicographic", "co-normal", "disjunctive", "OR",
"tensor", "Kronecker", "categorical", "direct"
)
}:= "Cartesian",
{
vertex_format::{string, symbol}:= "%A:%A",
labeltype::identical(string, symbol):= string
}
)
local
V:= GraphTheory:-Vertices~([G,H]),
E:= op~(4, [G,H]),
i, j, P:= [seq](seq([i,j], j= 1..nops(V[2])), i= 1..nops(V[1])),
Vgh:= (
curry(`if`(labeltype=string, sprintf, nprintf), vertex_format)
@ ((i,j)-> (V[1][i],V[2][j])) @ op
)~(P),
J:= op~(table(Vgh=~ P)),
Edge:= (u,v,k)-> J[v][k] in E[k][J[u][k]],
Eq:= (u,v,k)-> J[v][k] = J[u][k],
Rs:= table([
1= ((u,v)->     Edge(u,v,1) and       Eq(u,v,2)),
2= ((u,v)-> not Edge(u,v,1) and       Eq(u,v,2)),
3= ((u,v)->     Edge(u,v,1) and     Edge(u,v,2)),
4= ((u,v)-> not Edge(u,v,1) and     Edge(u,v,2)),
5= ((u,v)->     Edge(u,v,1) and not Edge(u,v,2)),
6= ((u,v)-> not Edge(u,v,1) and not Edge(u,v,2)),
7= ((u,v)->       Eq(u,v,1) and     Edge(u,v,2)),
8= ((u,v)->       Eq(u,v,1) and not Edge(u,v,2))
]),
Rels:= table([
("Cartesian", "box")=~ ``({1,7}),
("strong", "normal", "AND")=~ ``({1,3,7}),
"lexicographic"= ``({1,3,5,7}),
"modular"= ``({3,6}),
"homomorphic"= ``({5,7,8}),
("co-normal", "disjunctive", "OR")=~ ``({1,3,4,5,7}),
("tensor", "Kronecker", "categorical", "direct")=~ ``({3})
])
;
RelationGraph(
subs(_R= op(Rels[Rel]), (u,v)-> u<>v and ormap(k-> Rs[k](u,v), _R)),
Vgh
)
end proc
:
#Usage:
G:= GraphTheory:-CycleGraph(5): H:= GraphTheory:-CycleGraph(4):
GH:= GraphProducts(G, H, "co-normal", labeltype= symbol);
GH := Graph 3: an undirected unweighted graph with 20 vertices and 140 edge(s
```

The third argument can be any subset of {1,2,3,4,5,6,7,8}; seven of the subsets are given string names (with some synonyms) for convenience. The labeltype of the product's vetrtices defaults to string.

## ||...

What you want can be done with the || operator, which is very similar to cat, but without the evaluation issue that's causing your problem.

x||(a,b,c):= y =~ x^~(1,2,3);
rhs~([x||(a,c)])[];

## Can't be done...

Of course it can't be done. For it to be possible at all, at least the positions of the 0s must be the same.

 5 6 7 8 9 10 11 Last Page 7 of 361
﻿