## dharr

Dr. David Harrington

## 6416 Reputation

20 years, 23 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

## particular y(0)...

Yes, it seems the explicit solution has chosen y(0) = sqrt(2)/2 rather than an arbitrary value.

 > restart;
 > interface(version);

 > Physics:-Version();

 > libname;

 > restart;
 > ode:=y(x)*diff(y(x),x\$2)+diff(y(x),x)^2+1=0; IC:=D(y)(0)=1;

General solution has two constants as expected

 > sol:=dsolve(ode,explicit);

Requiring D(y)(0)=1 gives one equation that must be satisfied by c__1 and c__2. For the first solution this is

 > eval(diff(sol[1],x),x=0); eq1:=rhs(%)=1;

Let's use this to eliminate c2. We find a unique c__1 for each y(0) as expected.

 > isolate(eq1,c__2); s1:=eval(sol[1],%); eval(%,x=0);

If we take c__1 as -sqrt(2)/2 then the solution is

 > eval(s1,c__1=-sqrt(2)/2);

which is the same as the explicit solution

 > sol1:=dsolve([ode,IC],explicit);

 > odetest(sol1,[ode,IC,y(0)=sqrt(2)/2]);

But of course we could have taken any negative value for c__1, e.g. -1 (negative to make D(y)(0) positive, with positive c__2)

 > eval(s1, c__1=-1); odetest(%,[ode,IC,y(0)=1]);

Replace

`a[i, j] := -add(a[i, k], k = 1 .. N, k<>i) `

with

`a[i, j] := -(add(a[i, k], k = 1 .. N) - a[i, i])`

That is, just add all the terms and then remove the term you don't want.

## Degrees package or Units...

Here are two ways:

```with(Degrees):
evalf(17*sind(34)/sind(115));```

or

```with(Units):
evalf(17*sin(34*Unit(degree))/sin(115*Unit(degree)));```

(You can enter the units above with the Units palette.)

For degree minute seconds I think you will have to do the manipulations yourself

## eval...

As the error message says, the 2nd argument to eval should be an equation, not just the Matrix, i.e., eval(O, something = M)

## like this?...

Here the derivative of X__2 in terms of X__1 and its derivatives. I don't think it will ever be simple.

compact_derivatives.mw

## run worksheet...

If I run your worksheet, I get different answers than shown; u[3] etc now have their correct values.

As a separate issue, then you set a value for (u[i])[1], which doesn't make sense - now you make tables u[2] etc that are different from the table u, and you undo the assignments you made before.

## not autonomous...

"A system of two first order differential equations produces a direction field plot, provided the system is determined to be autonomous. In addition, a single first order differential equation produces a direction field (as it can always be mapped to a system of two first order autonomous differential equations). A system is determined to be autonomous when all terms and factors, other than the differential, are free of the independent variable. "

## no solutions in range...

The syntax is correct, just there are no solutions in this range. Works for -2..2

 >
 >

 >

 >

 >

## One way...

Seems like overkill, but here is one way. As for why, ...
(Note it is complexfreqvar that was s, not statevariable.)

 > restart;
 > foo:=proc()  local stemp;  if assigned(':-s') then    stemp := :-s;    :-s := ':-s';    DynamicSystems:-SystemOptions('complexfreqvar'=':-sv');    :-s := stemp;  else     DynamicSystems:-SystemOptions('complexfreqvar'=':-sv');  end if;  :-sv end proc:
 > s:="test"; foo(); s;

 >

## solve for z...

`sol := solve({eq1, eq2, eq3, x + y + z = 1}, {a, b, c, z}, explicit);`

## faster...

I have f(0, 12) taking 44 s, so you might get a result for f(0,30) in managable time.

 > restart

Use rationals

 > n1 := 423; x := 16; n2 := 81; y := 35; s1 := 1/10; b1 := 7; beta1 := 21/10; s2 := 1/10; b2 := 7; beta2 := 21/10;

 > A := (h, v, r) -> add(add(binomial(n1 - x, j)*binomial(n2 - y, l)*(-1)^(j + l)*GAMMA(s1 + h + 1)*GAMMA(s2 + v + 1)*(-1)^(h + v)*(beta1 - 1)^h*(beta2 - 1)^v/(h!*GAMMA(s1 + 1)*v!*GAMMA(s2 + 1)*(x + b1*s1 + j + b1*h + y + b2*s2 + l + b2*v)*(r + x + b1*s1 + j + b1*v)), j = 0 .. n1 - x), l = 0 .. n2 - y);

Warning, (in A) `l` is implicitly declared local

Warning, (in A) `j` is implicitly declared local

Result is a rational, but very small value, so will have numerical issues if you use floats

 > A(2,2,0);evalf(%);

Compare with a floating point calculation!

 >

 > B := (h, v, r) -> add(add(binomial(n1 - x, j)*binomial(n2 - y, l)*(-1)^(j + l)*GAMMA(s1 + v + 1)*GAMMA(s2 + h + 1)*(-1)^(h + v)*(beta1 - 1)^v*(beta2 - 1)^h/(v!*GAMMA(s1 + 1)*h!*GAMMA(s2 + 1)*(y + b2*s2 + l + b2*h + x + b1*s1 + j + b1*v)*(y + b2*s2 + l + b2*h - r)), j = 0 .. n1 - x), l = 0 .. n2 - y);

Warning, (in B) `l` is implicitly declared local

Warning, (in B) `j` is implicitly declared local

 > f := (r, upto) -> sum(sum(A(h, v, r) + B(h, v, r), v = 0 .. upto), h = 0 .. upto);

 > s := CodeTools:-Usage(f(0,12)):

memory used=14.57GiB, alloc change=214.03MiB, cpu time=32.14s, real time=44.15s, gc time=24.25s

 > s;

 > evalf(s);

 >

## agreed...

Maple likes the expanded form and one has to "fake it" in order to see exactly what you want. These sorts of manipulations can usually be done, but are not easy and IMO not worth the effort. Maple likes two fewer parentheses - isn't that good?

 >
 >

Use of %* instead of * stops the automatic multiplication

 >

Less ugly but more complicated

 >

 >

## indexing error...

Your code uses u in three different ways: as a procedure, as indexed variables, and as a 2-dimensional Array(1..40, 0..4). These should be using different symbols. By the time you enter eval_derivatives, u is the Array, and the other uses are gone. Then you use u[0], which means the zeroth row of the Array. But the rows start at 1, so you get the error message.

Also, I see later in eval_derivatives, you use T1i[i, k - j], when T1i is a 1D Array.

## generic...

By default, Maple assumes generic values of the coefficients. For special values, one needs to work harder. Here is one way to do it.

Edit: solve([eq1, eq2], [A__1, A__2, w], explicit) works also.

```SolveTools:-PolynomialSystem([eq1, eq2], [A__1, A__2, w], explicit)
```
 >
 > ;

If we treat it as a polynomial system

 >

 >

## table...

I rewrote this as a table. You don't say if you want all possible p and q sets for each expression, or just one. The following shows both possibilities. As I told you before, you were not properly detecting duplicates because the same expression was appearing in slightly different forms. evalc is enough to make a canonical form here; simplify is more obvious but too expensive.

There are various changes that could be make depending on how you want the output.

 > restart;
 > result:=table(sparse=NULL): #sparse=NULL only required if accumulating phi := (p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi) -> evalc((p__1*cosh(q__1*xi) - p__2*sinh(q__2*xi))/(p__3*cosh(q__3*xi) + p__4*sinh(q__4*xi))); for p__1 in [1, -1, I, -I] do     for p__2 in [1, -1, I, -I] do         for p__3 in [1, -1, I, -I] do             for p__4 in [ 0] do                 for q__1 in [ I, -I] do                     for q__2 in [ I, -I] do                         for q__3 in [1] do                             for q__4 in [ I, -I] do                                 # next line overwrites if we get the same result                                 result[phi(p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi)] := [p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4];                                 # OR, next line accumulates all possibilities                                 #result[phi(p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi)] ,= [p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4];                             od:                        od:                     od:                 od:             od:          od:      od: od:

All the expressions

 > exprns:=[indices(result,'nolist')]; print("Number of possibilities is ",numelems(exprns));

Find the corresponding p and q values. For just one

 > result[(-cos(xi) + sin(xi))*I/cosh(xi)]; result[exprns[3]];

For all of them

 > [entries(result,pairs)];

 >