dharr

Ask a Question

Create a Post

Name: Dr. David Harrington

22 Badges

Member for: 20 years, 335 days

Company/Institution: University of Victoria

Title: Professor or university staff

Location: Victoria, British Columbia, Canada

Website: http://web.uvic.ca/~dharr

Social Networks and Content at Maplesoft.com

Maple Application Center

Contact Dr. David Harrington

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.

MaplePrimes Activity

These are answers submitted by dharr

vertically aligned subscript and supersc...

Answered: dharr 8200

September 07 2024

1 1

Indexed subscripts to the power of a superscript are aligned vertically:

>	[f(x)][a]^b;

Download align.mw

convert to string...

Answered: dharr 8200

September 07 2024

0 0

Download Count_the_number_of_operations_in_A.mw

for loops or equivalent...

Answered: dharr 8200

September 07 2024

1 1

With for loops

n := 5;
A := Matrix(n, n);
for i to n do
    for j to n do
        A[i, j] := i*(1 + j);
    end do;
end do;
A;

Equivalent method

B := Matrix(n, n, (i, j) -> i*(1 + j));

Matrix.mw

(You should use a Matrix rather than an Array if you are going to be doing Matrix multiplication and other Matrix operations.)

multiplication...

Answered: dharr 8200

September 06 2024

0 0

Some of your multiplications are actually "." (Matrix multiplication, appears at bottom of line) instead of "*" (appears as centered dot).

Remove vs Statistics:-Remove...

Answered: dharr 8200

September 03 2024

3 0

The error message mentions Statistics:-Remove, which is the currently active Remove command after you loaded the Statistics package. To refer to the regular Remove command, you need to use the syntax

:-Remove(data, red)

save...

Answered: dharr 8200

September 02 2024

1 4

Use the save command to save d_remset to an .mpl file, which can be opened with a text editor, e.g.

save d_remset, cat(currentdir(),"/d_remset.mpl");

task sizes...

Answered: dharr 8200

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 8200

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 8200

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 8200

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 8200

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 8200

August 12 2024

1 1

Here are two ways:

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

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 8200

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 8200

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 8200

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

First