dharr

task sizes...

Answered: dharr 8440

September 01 2024

0 3

Based on @Carl Love's, comment, I rewrote makeproc to not access the global Equation and j11, and also to not use Vars and Vars1 since I was worried about memory. It turns out that does help go to N=160 and is also faster, but the main issue actually seems to be the interaction between the value of N and the "task size" number, which was 200 in @vsubramanian's code. As you increase N it seems you have to increase this number. I'm not sure what his means; maybe you spawn too many tasks or maybe it is a memory issue. For N=160, I used 4000 as in the code below. But this wasn't always consistent, perhaps because of some memory issues, and might be easier if I wasn't comparing the old and new code. And it means the benefits of multthreading are not so great, at least on this machine with numcpus = 12.

I hacked quite a bit with the makeproc code, but I think it is functionally the same as before. It seems the variable numbers in the first part (making Equation) and second part (differentiation) are the same, but in case that was specific to this example I didn't mess with it, as per @vsubramanian's PS.

>

restart:

>	Digits:=15;

Original

>

makeproc:=proc(n0,nf,Nt,Eqode::Array,Vars::list,Vars1::list,Equation::Array,j11::Matrix(storage=sparse))
local i,j,LL,LL2,L,i1,varsdir,varsdir1,node,eqs;

for i from n0 to nf do
eqs:=rhs(Eqode[i]):
L:=indets(eqs)minus{t};
if nops(L)>0 then
LL := [seq(ListTools:-Search(L[j], Vars),
             j = 1 .. nops(L))];
LL2 := ListTools:-MakeUnique(LL);
if LL2[nops(LL2)] = 0 then
             LL := [seq(LL2[i1], i1 = 1 .. nops(LL2) - 1)];
         end if;
         if LL2[1]= 0 then LL:=[seq(LL2[i1],i1=2..nops(LL2))]: end:#print(1,LL);
varsdir:=[seq(Vars[LL[i1]]=0.5*uu[LL[i1]]+uu_old[LL[i1]],i1=1..nops(LL))];
else varsdir:=[1=1]:end:

Equation[i]:=uu[i]-deltA*subs(varsdir,t=tnew,eqs):#print(i,Equation[i]);
        L := indets(Equation[i]);
        LL := [seq(ListTools:-Search(L[j], Vars1),j = 1 .. nops(L))];
        LL2 := ListTools:-MakeUnique(LL);
        if LL2[nops(LL2)] = 0 then
            LL := [seq(LL2[i1], i1 = 1 .. nops(LL2) - 1)];
        end if;
        if LL2[1]= 0 then LL:=[seq(LL2[i1],i1=2..nops(LL2))]: end:

        for j to nops(LL) do
            j11[i, LL[j]] := diff(Equation[i], uu[LL[j]]);
        end do;

od:
end proc:

New version assumes the variable ordering in Vars and Vars1. Output to tables (not global array/matrix) from which matrices can be constructed

>

makeproc3:=proc(n0,nf,N,Eqode::Array)
local i, L, varsdir, eqs, Equation, j11, varnum, subst, uvar, cvar;
varnum := (i,j)->(j-1)*N+i;
subst := proc(var) local vn; # takes var = c[i,j](t) and makes c[i,j](t)=0.5*uu[k]+uu_old[k]
vn := varnum(op(op(0,var))); # get variable number from subscripts [i,j]
var = 0.5*'uu'[vn]+'uu_old'[vn];
end proc;

Equation := table();
j11 := table();

for i from n0 to nf do
eqs := rhs(Eqode[i]):
L := indets(eqs) minus {'t'};

varsdir := [seq(subst(cvar), cvar in L)];

Equation[i] := 'uu'[i]-'deltA'*subs(varsdir,'t'='tnew',eqs):
        L := indets(Equation[i],'specindex'('uu'));

        for uvar in L do
            j11[i, op(uvar)] := diff(Equation[i], uvar);
        end do;

end do;
[eval(Equation),eval(j11)]
end proc:

>

>	CodeTools[ThreadSafetyCheck](makeproc); CodeTools[ThreadSafetyCheck](makeproc3);

>	#infolevel[all]:=10;

>

>	N:=160;h:=1.0/N:

>	with(Threads[Task]):

>	Eqs:=Array(1..N,1..N):

>

for i from 1 to N do for j from 1 to N do
if i = 1 then left:=0: else left:=(c[i,j](t)-c[i-1,j](t))/h:end:
if i = N then right:=-0.1: else right:=(c[i+1,j](t)-c[i,j](t))/h:end:
if j = 1 then bot:=0: else bot:=(c[i,j](t)-c[i,j-1](t))/h:end:
if j = N then top:=-0.1: else top:=(c[i,j+1](t)-c[i,j](t))/h:end:
Eqs[i,j]:=diff(c[i,j](t),t)=(right-left)/h+(top-bot)/h-c[i,j](t)^2:
od:od:

>

>	eqs1:=Array([seq(seq(Eqs[i,j],i=1..N),j=1..N)]):

>	ics:=[seq(seq(c[i,j](t)=1.0,i=1..N),j=1..N)]:

>	Vars:=[seq(seq(c[i,j](t),i=1..N),j=1..N)]:

>	Vars1:=[seq(uu[i],i=1..N^2)]:

>	Equation:=Array(1..N^2):j11:=Matrix(1..N^2,1..N^2,storage=sparse):eqs:=copy(Equation):

>	CodeTools:-Usage(makeproc(1,N^2,N^2,eqs1,Vars,Vars1,Equation,j11));

memory used=375.28MiB, alloc change=56.00MiB, cpu time=3.98s, real time=6.77s, gc time=1.20s

>

j11[5,5];

Save original results for comparison

>	j11copy:=copy(j11): Equationcopy:=copy(Equation):

Run new routine

>	tbls:=CodeTools:-Usage(makeproc3(1,N^2,N,eqs1)):

memory used=175.28MiB, alloc change=-16.00MiB, cpu time=2.11s, real time=2.52s, gc time=875.00ms

Construct Array and Matrix from the tables

>	Equation := Array(1..N^2, tbls[1]): j11 := Matrix(1..N^2,1..N^2, tbls[2], storage = sparse):

Check results are the same

>	j11[5,5]; EqualEntries(Equation, Equationcopy); EqualEntries(j11, j11copy);

Continuation routine merges the tables - tablemerge not listed as thread safe

>	cont := proc(tbls1, tbls2) # merge the tables [table([entries(tbls1[1],'pairs'),entries(tbls2[1],'pairs')]), table([entries(tbls1[2],'pairs'),entries(tbls2[2],'pairs')])]; end proc:

>	CodeTools[ThreadSafetyCheck](cont);

>

makeprocDistribute := proc(i_low, i_high, N, Eqode::Array)
local i_mid;
if 4000 < i_high - i_low then
i_mid:=iquo(i_low+i_high,2):
Threads:-Task:-Continue(cont,
Task = [makeprocDistribute, i_low, i_mid, N, Eqode],
Task = [makeprocDistribute, i_mid + 1, i_high, N, Eqode]);
else
return makeproc3(i_low, i_high,N,Eqode);
end if;
end proc:

>

j11[5,5];

>

N^2;

>

NN:=N^2;

>	tbls:=CodeTools:-Usage(Start(makeprocDistribute,1,NN,N,eqs1)):

memory used=214.42MiB, alloc change=347.19MiB, cpu time=10.89s, real time=2.12s, gc time=20.06s

Construct Array and Matrix from the tables

>	Equation := Array(1..N^2, tbls[1]): j11 := Matrix(1..N^2,1..N^2, tbls[2], storage = sparse):

Check results are the same

>	j11[5,5]; EqualEntries(Equation, Equationcopy); EqualEntries(j11, j11copy);

>

Download makeproctest3.mw

Edit: new code only, success with N=500 (system is faster b/c on high performance setting):
makeproctestnew.mw

min and max...

Answered: dharr 8440

August 30 2024

2 6

KR is a table, so you can use min(entries(KR),'nolist'). But KR[N], h0 and tau0 (not tao0?) are just values, so it's not clear what you mean.

In general if you are calculating a value in a loop you can do something like to to calculate the minimum:

minval:=infinity;
for i to 100 do
  val := ...; # some calculation
  if val < minval then minval := val end if;
end do;

and similarly for the maximum.

Or you could collect all the values in a Vector V and just use min(V); (or use a list or Array).

odeplot code...

Answered: dharr 8440

August 30 2024

1 3

odeplot handles the different types of output, so you could look at the code there as a starting point - see lines 5-14 in showstat(`plots/odeplot`).

particular y(0)...

Answered: dharr 8440

August 26 2024

1 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;

$"C:\Program Files\Maple 2024\lib"$

>

restart;

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

y(x)*(diff(diff(y(x), x), x))+(diff(y(x), x))^2+1 = 0

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;

$eval(diff(y(x), x), {x = 0}) = -(1/2)*2^(1/2)*c__1/c__2^(1/2)$

-(1/2)*2^(1/2)*c__1/c__2^(1/2) = 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);

c__2 = (1/2)*c__1^2

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]);

Download dsolve_constants.mw

add and then subtract...

Answered: dharr 8440

August 23 2024

2 0

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

Answered: dharr 8440

August 12 2024

1 1

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

Answered: dharr 8440

August 11 2024

0 0

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

Answered: dharr 8440

August 09 2024

2 6

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

Answered: dharr 8440

August 09 2024

1 1

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.

Download table.mw

not autonomous...

Answered: dharr 8440

August 08 2024

3 1

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

Answered: dharr 8440

August 08 2024

1 0

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

Download fsolve.mw

One way...

Answered: dharr 8440

August 08 2024

1 0

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

Download dynamic.mw

solve for z...

Answered: dharr 8440

August 08 2024

1 1

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

faster...

Answered: dharr 8440

August 08 2024

1 2

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

For future reference please upload your worksheet using the green up-arrow in the Mapleprimes editor.

>

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

proc (h, v, r) local l, j; options operator, arrow; add(add(binomial(n1-x, j)*binomial(n2-y, l)*(-1)^(l+j)*GAMMA(s1+h+1)*GAMMA(s2+v+1)*(-1)^(h+v)*(beta1-1)^h*(beta2-1)^v/(factorial(h)*GAMMA(s1+1)*factorial(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) end proc

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

proc (h, v, r) local l, j; options operator, arrow; add(add(binomial(n1-x, j)*binomial(n2-y, l)*(-1)^(l+j)*GAMMA(s1+v+1)*GAMMA(s2+h+1)*(-1)^(h+v)*(beta1-1)^v*(beta2-1)^h/(factorial(v)*GAMMA(s1+1)*factorial(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) end proc

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

proc (r, upto) options operator, arrow; sum(sum(A(h, v, r)+B(h, v, r), v = 0 .. upto), h = 0 .. upto) end proc

>	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);

>

Download sum.mw

agreed...

Answered: dharr 8440

August 07 2024

2 1

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?