## dharr

Dr. David Harrington

## 4791 Reputation

19 years, 34 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

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

## and the loop is not right...

aside from the main issue raised by @rlopez your loop should be "from 0 to 99" not "by 0 to 99"

## fsolve works for me...

fsolve(y, x = 0 .. 2e6);

works for me.

## A special case?...

@leslieglasser Still not clear what the possible inputs are. If it is a monic polynomial and the other coefficents are the variables (as in your example), then the following works:

 > restart;with(LinearAlgebra):
 (1)
 > M1:=:M2:=^(%T):
 > M:=<,M2>; Determinant(M);
 (2)
 >

I don't have any suggestions for more complicated cases.

@leslieglasser  Well you can make a 3x3 one by bordering the 2x2 one with zeroes except for one in the 1,1 entry. But I'm guessing that isn't what you want. Can you be more specifc about your exact requirements? - I don't know what separated variables mean or nomograph in this context.

## closer?...

int(Heaviside(1-x^2-y^2), [x=0..1,y=0..1]); #should be Pi/4

gives -Pi/4

(I thought using even symmetry might help, but it went negative.)

## Set x...

@dharr If you set x=Pi/2 rather than solving for it, then you get the a, b, c values that you want.

## Not correct solution...

ab should be a*b, but if I've interpreted your equation correctly, then odetest suggests that your solution is not correct.

diff(psi(x), x, x)-((b^2-a*(2*p+3))*x^2+2*a*b*x^4+a^2*x^6-energy)*psi(x) = 0;
schro := eval(%,energy=(2*p+1)*b);
odetest(psi(x)=(x^p)*exp((-a*(x^4))/4)-(b*(x^2)/2),schro);

does not return zero, meaning it is not a solution. Using dsolve gives something with Heun functions, but you probably need to check you have the right equation first.

## complicated...

@sunit I'm not entirely sure what you mean. You could certainly write a procedure that tested whether there were minus signs to its argument and adjusted things before passing to sin. But covering all the unexpected cases (the above one was certainly unexpected for me; I just mistakenly typed Phi instead of phi) is probably going to take more effort than it is worth. As I noted, Maple isn't human, and often outputs things in forms that are unexpected (but mathematically correct)

## Maple has a mind of its own...

 > restart;
 > q1:=sin(omega*(T0-T)+Phi);
 (1)
 > map(expand,q1);
 (2)
 > q2:=sin(omega*(T0-T)+phi);
 (3)

Notice the Maple has changed rearranged to have a minus before the sin.

 > map(expand,q2);
 (4)

map applies expand to each operand. For q1, the only operand is the argument to sin, but Maple considers q2 as a product of two operands, the -1 and the sin, so it expands both

 > whattype(q1);op(q1);
 (5)
 > whattype(q2);op(q2);
 (6)

To get what you want requires you to operate the expand on the (first) operand of the second operand, which can be done as follows

 > applyop(expand,[2,1],q2);
 (7)

But this is a rather complicated solution, which requires attention in individual cases. Of course you could just stick with Phi. In general, getting Maple to rearrange things into the exact form you want is difficult. Depending on your ultimate objective, it may not be worth it. But perhaps others have a simpler solution.

 >

## OK...

I misread the correct answer at line (13) followed by the comment "The results are worse than single iteration using jacobian" to mean that you didn't think the default solver was working. I see you have now clarified that with an edit.

As for correctly coding lsode I don't have much comment; I have tried to move away low level coding.

## Maple doesn't know this...

FunctionAdvisor(hypergeom([a1,a2],[b1,b1],z)); returns "The system is unable to compute the "aymptotic expansion" for hypergeom([a1,a2],[b1,b2],z)"

The  NIST handbook gives a formula but it's pretty messy:

http://dlmf.nist.gov/16.11#E7

## Rosenbrock seems to be working...

I programmed mgear in FORTRAN long ago, and then used it in Maple, then Maple's lsode, then Maple's rosenbrock, and I haven't noticed any difficulties.

I can't tell exactly what you are doing, but it seems as if you know what the solution should be? What is that solution and why do you think the 0.3496  answer is wrong? There doesn't seem to be anything pathological about the DE. At least locally, the answer is right, as this shows:

 > restart;
 > de:=diff(y(x),x)=-y(x)^2*exp(-y(x))-10*y(x)*(1+0.01*exp(y(x)));
 (1)

Default (non-stiff) solver)

 > ans:=dsolve({de,y(0)=1},y(x),numeric); ans(0.1);
 (2)

Lsode solver

 > ans:=dsolve({de,y(0)=1},y(x),numeric,method=lsode); ans(0.1);
 (3)

Default stiff solve (Rosenbrock)

 > ans:=dsolve({de,y(0)=1},y(x),numeric,stiff=true); ans(0.1);
 (4)

Check out the answer from Rosenbrock, which is easier to do with the listprocedure output. The numerical derivative at 0.1 is equal to the value of the right-hand side of the de to four decimal places

 > ans:=dsolve({de,y(0)=1},y(x),numeric,stiff=true,output=listprocedure); yx:=rhs(ans[2]); yx(0.1); fdiff(yx(x),x=0.1); eval(rhs(de),y(x)=yx(0.1));
 (5)
 >

## Worksheet and choice of 1-D or 2-D...

I teach a short introduction to Maple course for graduate students in chemistry, and had to think carefully about this. I decided that since the point of the course is to teach them how to do math and simple programming with Maple rather than optimize the way the document looks, the simple logic of the worksheet mode was more suitable than document mode. 1-D or 2-D was a harder choice, and I let them submit in either. Some like 2-D and some 1-D. Those who choose 2-D seem to stick with it, even though it trips them up on occasion.

On the other hand, I demonstrate it all in 1-D. They will be looking at 1-D when displaying procedures, and if they end up doing any significant programming they will likely use it. I suggest that even if using 2-D they may wish to enter multiplication as * to avoid the errors of missing spaces, and to learn the 1-D notation at the same time.

I don't recommend extended typesetting, since seeing output with BesselJ will remind them how to enter that Bessel function. Like document mode, it is a feature for more advanced users, that they can teach themselves.

@Christopher2222 So if you are going to use the geometry package a lot, you can load all the commands using

with(geometry);

## from bits to bits...

@adrian5213 Not sure what you mean by a "binary number" (a string of 0's and 1's? a sequence? a list?). If you are OK with a list of 0's and 1's then the following works.

n:=[1,0,1,0,1,1,0,0,1];
ListTools:-Reverse(n): #reverse order since Bits package works in lsb first order
Bits:-Join(%): #convert to internal (integer) format
[Bits:-GetBits(-%,-1..0,bits=8)]: #take the negative and convert back to a list
subsop(1=0,%); #force the first bit to zero

produces

[0, 0, 1, 0, 0, 1, 1, 1]

I assumed the input had first bit 1, but you would probably want to check.

This of course gets Maple to do the work. If you want to do it "by hand", then you could just construct the usual algorithm with list entries.

 First 38 39 40 41 42 43 44 Page 40 of 45
﻿