dharr

Dr. David Harrington

8270 Reputation

22 Badges

20 years, 356 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
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

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

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

restart;

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 11, June 25 2024 Build ID 1835466`

Physics:-Version();

`The "Physics Updates" package is not installed`

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

(D(y))(0) = 1

General solution has two constants as expected

sol:=dsolve(ode,explicit);

y(x) = (-2*c__1*x-x^2+2*c__2)^(1/2), y(x) = -(-2*c__1*x-x^2+2*c__2)^(1/2)

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

y(x) = (c__1^2-2*c__1*x-x^2)^(1/2)

y(0) = (c__1^2)^(1/2)

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

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

y(x) = (1/2+2^(1/2)*x-x^2)^(1/2)

which is the same as the explicit solution

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

y(x) = (1/2+2^(1/2)*x-x^2)^(1/2)

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

[0, 0, 0]

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

y(x) = (-x^2+2*x+1)^(1/2)

[0, 0, 0]

 

Download dsolve_constants.mw

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.

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

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

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

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

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

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

restart

f := sin(x[1]+x[2])-exp(x[1])*x[2]; g := x[1]^2-x[2]-2; cond := {seq(x[i] = -1 .. 1, i = 1 .. 2)}; fsolve({f, g}, cond)

sin(x[1]+x[2])-exp(x[1])*x[2]

x[1]^2-x[2]-2

{x[1] = -1 .. 1, x[2] = -1 .. 1}

fsolve({sin(x[1]+x[2])-exp(x[1])*x[2], x[1]^2-x[2]-2}, {x[1], x[2]}, {x[1] = -1 .. 1, x[2] = -1 .. 1})

plot3d([f, g], x[1] = -2 .. 2, x[2] = -2 .. 2)

cond2 := {seq(x[i] = -2 .. 2, i = 1 .. 2)}; fsolve({f, g}, cond2)

{x[1] = -2 .. 2, x[2] = -2 .. 2}

{x[1] = -.6687012050, x[2] = -1.552838698}

NULL

Download fsolve.mw

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;

"test"

sv

"test"

 

Download dynamic.mw

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

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;

423

16

81

35

1/10

7

21/10

1/10

7

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

11727199469241514779864393729287999595266717038898384405454822453889989726964576758757729551099323484694782170822448350244634045448648056150840927233527350118380124119264858615283201523602923810903830006171445640638533014964609388484372553536841310967853824558690932280088151544584971199415788181111551306532754880119508566210291084045598562471540440317095151211911113829813243986384020095260781211101363284506805774438303489939264952762167438202379620004245323297602579179457497328410039101743209644168235298998961095392278935981336337925654946855818892100272459692432824817768478928223155226897063053135874218725475831316965225841532156062303907670391646481346412180543470550738955412983380825233032396184376956505282354537912590976013025195012646732190615857238658037431378033943474292755126953125/6382265268317256219633710058998499818340873590563036702803822950484199699636302738462052037083809489945501654638778705432982835407262106422716055822223400225220629848888603203921980303115472012139298206087469920885254016501522878034151924440738350864243996392880003254552368644323922787997879524243559004914030851003740366052093558872264913393738584398643237741388376227019971696955206515106578177861348809414199987981935214749218771677333981945489142846776730080142869549055976924435701756372241817375336920507491648219754129346141714754707551648893025832772408084923821700316568058764587634659823042662378587917652228685077432247136823506786967161644710468316981814899035776299987989971099597861998466613316946045883428522261542320902882571202585668513207901611181230727412654811333329187224660577277205763279177427978308304241574151519438327049907892296647014672715620294656

0.1837466632e-77

Compare with a floating point calculation!

A(2, 2, 0.)

-0.3955624180e121

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;

114249022902597666588491400529505235216742049616231180869631619194167698464302221521109360278934161699249264126826236228253701293270492644550135177145016349814982113455183820820735876815796463679992519208283503642450130209375877055834771405738601458875577451563602032933837943025767731232207854934759236939341322568521090401156938484109405395425106996068885180201446511126644125409694777192264543924174126949778005106553840262917702653247531868774015140088404720078587672388076732410686714770939922375245824000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000/4478112673625161634402912271274057087521425782459964471363322850631457597742398047074088546911131451856756368982860023346654177624065136219112870362974132591633772922976976016146648645581986021657025573763585023497503972196876316120665545360282282680322522472194461718178716824852465782992306776138047404097571796261065850414338472933710797700189622470162561730240651195333635833352277583626670180879730361251512986788399320104862652482661166302364468264845833967133426250993316566352774199352462247496894848151055424852777184106477710417868463737767180212637570654326336497085984319988974445914462490654328034321182352155451367268612522705340661230039301539950445376032066380305314154637436491880112866704742279027028488977660000585169277423676424285496535951865146024850581707926633898571343437487693846692653268159361127272552016275556640003360341921899571661034103268149587939236207824679329125781964757872174171342424591595209864084608515982176191174973912465653082613424079377178453298484585210385316658608895204147668705453968995722653400843313999068483592478297363991

evalf(s);

0.2551276201e-55

NULL

Download sum.mw

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?

restart

p := 9*csc(theta)^2+9

9*csc(theta)^2+9

Use of %* instead of * stops the automatic multiplication

q := content(p)*%*primpart(p); value(q)

`%*`(9, csc(theta)^2+1)

9*csc(theta)^2+9

Less ugly but more complicated

InertForm:-Display(content(p)*%*primpart(p), inert = false)

0, "%1 is not a command in the %2 package", _Hold, Typesetting

NULL

Download content.mw

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.

First 14 15 16 17 18 19 20 Last Page 16 of 82