Advanced Search

dharr

Dr. David Harrington

8270 Reputation

22 Badges

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

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

LinearSolve...

dharr 8270
July 28 2022
1 3

Since the system is linear for the variables you want to solve for (for a given t), you can use LinearSolve.

> 

restart
> 

with(LinearAlgebra):
> 

EQ := Matrix(4, 1, {(1, 1) = 32.1640740637930*Tau[1]-0.172224519601111e-4*Tau[2]-0.270626540730518e-3*Tau[3]+0.1570620334e-9*P[1]+0.3715450960e-14*sin(t), (2, 1) = -0.172224519601111e-4*Tau[1]+32.1667045885952*Tau[2]+0.587369829416537e-4*Tau[3]-0.1589565489e-8*P[1]+0.1004220091e-12*sin(t), (3, 1) = -0.270626540730518e-3*Tau[1]+0.587369829416537e-4*Tau[2]+32.1816411689934*Tau[3]-0.7419658527e-8*P[1]+0.5201228088e-12*sin(t), (4, 1) = 0.1570620334e-9*Tau[1]-0.1589565489e-8*Tau[2]-0.7419658527e-8*Tau[3]+601.876235436204*P[1]})
> 

EQlist := convert(EQ, list);

[32.1640740637930*Tau[1]-0.172224519601111e-4*Tau[2]-0.270626540730518e-3*Tau[3]+0.1570620334e-9*P[1]+0.3715450960e-14*sin(t), -0.172224519601111e-4*Tau[1]+32.1667045885952*Tau[2]+0.587369829416537e-4*Tau[3]-0.1589565489e-8*P[1]+0.1004220091e-12*sin(t), -0.270626540730518e-3*Tau[1]+0.587369829416537e-4*Tau[2]+32.1816411689934*Tau[3]-0.7419658527e-8*P[1]+0.5201228088e-12*sin(t), 0.1570620334e-9*Tau[1]-0.1589565489e-8*Tau[2]-0.7419658527e-8*Tau[3]+601.876235436204*P[1]]

> 

V := convert(Matrix(1, 4, {(1, 1) = Tau[1], (1, 2) = Tau[2], (1, 3) = Tau[3], (1, 4) = P[1]}), list)

[Tau[1], Tau[2], Tau[3], P[1]]

> 

A, b := GenerateMatrix(EQlist, V)

A, b := Matrix(4, 4, {(1, 1) = 32.1640740637930, (1, 2) = -0.172224519601111e-4, (1, 3) = -0.270626540730518e-3, (1, 4) = 0.1570620334e-9, (2, 1) = -0.172224519601111e-4, (2, 2) = 32.1667045885952, (2, 3) = 0.587369829416537e-4, (2, 4) = -0.1589565489e-8, (3, 1) = -0.270626540730518e-3, (3, 2) = 0.587369829416537e-4, (3, 3) = 32.1816411689934, (3, 4) = -0.7419658527e-8, (4, 1) = 0.1570620334e-9, (4, 2) = -0.1589565489e-8, (4, 3) = -0.7419658527e-8, (4, 4) = 601.876235436204}), Vector(4, {(1) = -0.3715450960e-14*sin(t), (2) = -0.1004220091e-12*sin(t), (3) = -0.5201228088e-12*sin(t), (4) = 0})

> 

soln := Equate(V, LinearSolve(A, b));

[Tau[1] = -0.1156532164e-15*sin(t), Tau[2] = -0.3121894613e-14*sin(t), Tau[3] = -0.1616209236e-13*sin(t), P[1] = -0.2074537757e-24*sin(t)]

For t = 1

> 

eval(soln, t = 1.);

[Tau[1] = -0.9731882590e-16, Tau[2] = -0.2626983734e-14, Tau[3] = -0.1359993177e-13, P[1] = -0.1745663329e-24]

NULL

Download SoalNewton.mw

pointplot?...

dharr 8270
July 26 2022
0 0

It is not very clear what information you want to show, or what form you have the data in. There are two sets of data in the table, but only one set of heights in the image, and those heights are larger than the values in the table. The lines in between the points do not seem to mean anything. The grid of x1, x2 values is not complete - are heights for missing values to be set to zero?

If only the heights in one column, say the first one, are of interest, then you could use plot with style=point (or pointplot) to plot the heights at the various x1,x2 points.

Download plot.mw

solve...

dharr 8270
July 26 2022
1 1

To find which value(s) of f(x) make the equation a*f(x)=0 the same on both sides, use solve:
solve(a*f(x) = 0, f(x)) assuming a>0

gives 0 for f(x). Simplify will simplify each side of the equation separately, but does not know that you want to find f(x).

quo and rem...

dharr 8270
July 22 2022
2 1


 

> 

p:=x^3-3*x^2-12*x+19;
d:=x^2+x-6;

x^3-3*x^2-12*x+19

x^2+x-6

Quotient

> 

q:=quo(p,d,x,'r');

x-4

Remainder

> 

r;

-2*x-5

> 

expand(d*q+r);

x^3-3*x^2-12*x+19

> 

``


 

Download quo.mw

via Eigenvalues...

dharr 8270
July 15 2022
3 0
> 

restart
> 

with(LinearAlgebra):
> 

A := Matrix([[w-k^2-2*sqrt(2)*sqrt(b^2+d^2)*k+2*(b^2+d^2), 2*(b^2-d^2), 0, 2*b*d, 0, -2*b*d], [2*(b^2-d^2), -w-k^2+2*sqrt(2)*sqrt(b^2+d^2)*k+2*(b^2+d^2), 2*b*d, 0, -2*b*d, 0], [0, 4*b*d, w-k^2-2*sqrt(2)*sqrt(b^2+d^2)*k+2*(b^2+d^2), 2*d^2, 0, 2*b^2], [4*b*d, 0, 2*d^2, -w-k^2+2*sqrt(2)*sqrt(b^2+d^2)*k+2*(b^2+d^2), 2*b^2, 0], [0, -4*b*d, 0, 2*b^2, w-k^2-2*sqrt(2)*sqrt(b^2+d^2)*k+2*(b^2+d^2), 2*d^2], [-4*b*d, 0, 2*b^2, 0, 2*d^2, -w-k^2+2*sqrt(2)*sqrt(b^2+d^2)*k+2*(b^2+d^2)]]):
> 

evs := Eigenvalues(A)

evs := Vector(6, {(1) = 2*b^2+2*d^2-k^2+sqrt(-4*sqrt(2)*sqrt(b^2+d^2)*k*w+4*b^4+8*b^2*d^2+8*b^2*k^2+4*d^4+8*d^2*k^2+w^2), (2) = 2*b^2+2*d^2-k^2-sqrt(-4*sqrt(2)*sqrt(b^2+d^2)*k*w+4*b^4+8*b^2*d^2+8*b^2*k^2+4*d^4+8*d^2*k^2+w^2), (3) = 2*b^2+2*d^2-k^2+sqrt(-4*sqrt(2)*sqrt(b^2+d^2)*k*w+4*b^4+8*b^2*d^2+8*b^2*k^2+4*d^4+8*d^2*k^2+w^2), (4) = 2*b^2+2*d^2-k^2-sqrt(-4*sqrt(2)*sqrt(b^2+d^2)*k*w+4*b^4+8*b^2*d^2+8*b^2*k^2+4*d^4+8*d^2*k^2+w^2), (5) = 2*b^2+2*d^2-k^2+sqrt(-4*sqrt(2)*sqrt(b^2+d^2)*k*w+4*b^4+8*b^2*d^2+8*b^2*k^2+4*d^4+8*d^2*k^2+w^2), (6) = 2*b^2+2*d^2-k^2-sqrt(-4*sqrt(2)*sqrt(b^2+d^2)*k*w+4*b^4+8*b^2*d^2+8*b^2*k^2+4*d^4+8*d^2*k^2+w^2)})

There are two values of w that make an eigenvalue zero, and therefore the determinant zero:

> 

{seq(solve(evs[i], w), i = 1 .. 6)};

{(2*2^(1/2)*(b^2+d^2)^(1/2)-(-4*b^2-4*d^2+k^2)^(1/2))*k, (2*2^(1/2)*(b^2+d^2)^(1/2)+(-4*b^2-4*d^2+k^2)^(1/2))*k}

> 

simplify(Determinant(eval(A, w = w1)));

0

> 

simplify(LinearAlgebra:-Determinant(eval(A, w = w2)));

0

> 

NULL

NULL

Download determinant.mw

LinearAlgebra...

dharr 8270
July 06 2022
2 1

The Column command is in the LinearAlgebra package, so I'm guessing there wasn't with(LinearAlgebra) before it to load the package.

floor...

dharr 8270
July 05 2022
3 2
> 

plot(floor(x), x = -2 .. 5, discont = true)

> 

NULL

Download floor.mw

don't need to pair...

dharr 8270
June 21 2022
3 2

If you use plot with style=point, then you can use the lists directly without pairing.

plot(P,Q,style=point,symbol=solidcircle,symbolsize=15,color=red,view=[-1..1,0..30]);

 

Inert in Maple 2015...

dharr 8270
June 18 2022
1 1

In Maple 2015, the following works. (Not displaying correctly in Mapleprimes - see worksheet).
 

> 

restart;
> 

with(InertForm):
> 

a:=1: b:= 2: c:=1: d:= 6: e:= 2:
> 

P:=`%*`(a,b,c,d,e);

`%*`(1, 2, 1, 6, 2)

> 

Display(P);

`%*`(1, 2, 1, 6, 2)

> 

Display(P,inert=false);

1.2.1.6.2

> 

``

 

Download Inert.mw

nested coeff...

dharr 8270
June 18 2022
3 1


 

> 

eq:=2*cosh(xi)^8*alpha*B[2]^2 + 4*A[2]*sinh(xi)^6*cosh(xi)^2 + 4*A[1]*sinh(xi)^5*cosh(xi)^3 + 4*sinh(xi)^4*cosh(xi)^4*A[0];

2*cosh(xi)^8*alpha*B[2]^2+4*A[2]*sinh(xi)^6*cosh(xi)^2+4*A[1]*sinh(xi)^5*cosh(xi)^3+4*sinh(xi)^4*cosh(xi)^4*A[0]

> 

coeff(coeff(eq,cosh(xi)^2),sinh(xi)^6);

4*A[2]

 

Download coeff.mw

Matrix row by row...

dharr 8270
June 17 2022
3 0

This may not be the most efficient way to do it; a lot depends on if you have simple ways to generate the arguments or the functions or function names.
 

> 

restart;
> 

f1:=x->x^2;f2:=x->x^3;f3:=x->sin(x);

proc (x) options operator, arrow; x^2 end proc

proc (x) options operator, arrow; x^3 end proc

proc (x) options operator, arrow; sin(x) end proc

> 

nrows:=4:
A:=Matrix(nrows,3):
for i to nrows do
  A[i,...]:=<f1(i/2.)|f2(i*2.)|f3(i*1.5)>;
end do:
> 

A;

Matrix([[.2500000000, 8., .9974949866], [1.000000000, 64., .1411200081], [2.250000000, 216., -.9775301177], [4.000000000, 512., -.2794154982]])

> 

ExcelTools:-Export(A,cat(currentdir(),"/test.xlsx"));
> 

 

 

Download Excel.mw

one condition...

dharr 8270
June 17 2022
4 1

For a first order ode, you can specify only one initial/boundary condition. Then your first case returns a solution in terms of Heaviside.

In my version of Maple, I don't get a solution for the second one if I just take the phi(0) condition.

> 

restart;
> 

eq:=diff(theta(t),t)=-(A[1]*Dirac(t-T[1])+A[2]*Dirac(t-T[2]));

diff(theta(t), t) = -A[1]*Dirac(t-T[1])-A[2]*Dirac(t-T[2])

> 

bcon := theta(0)=-v[1];

theta(0) = -v[1]

> 

ans:=dsolve({eq,bcon},theta(t));

theta(t) = -A[1]*Heaviside(t-T[1])-A[2]*Heaviside(t-T[2])-v[1]+A[1]*Heaviside(-T[1])+A[2]*Heaviside(-T[2])

> 

 


 

Download Dirac.mw

scanf...

dharr 8270
June 16 2022
1 0

Not sure if you already have this in a string or that is the file contents. Either way, this should work.

> 

restart;
> 

str:="[[1. 0.\n]\n[0. 1.]]"

"[[1. 0. ] [0. 1.]]"

Scan to give a Matrix of size 2,2, ignoring characters "[", "]" and ".".

You could read this in directly from a file using fscanf instead of sscanf.

See ?rtable_scanf for more information

> 

A:=sscanf(str,"%{2,2;d([].)}am")[];

A := Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1})

> 

 

Download MatrixInput.mw

with finite initial velocity...

dharr 8270
June 15 2022
4 1

I didn't submit this earlier since I had the same problem that @Rouben Rostamian had. But if you assume that it doesn't start at rest, but is given an initial finite y-component of the velocity, then you can derive the semicubical parabola.
 

> 

restart

Solve isochrone problem (solution is Neile's semicubical parabola).

We start a mass m at positive coordinates (x__0, y__0) with x component of velocity zero and y component of velocity -v__y and let it move toward the origin under the action of gravity along some curve for which the y (height) coordinate changes by equal amounts in equal times. What is the curve?

 

Solution. Run the problem backwards. At time zero, let the mass leave the origin with velocity components v__x(0) = v__0 and v__y. The problem statement means that v__y is constant.

Let t be the time the particle arrives at coordinates (x,y), so we know y = v__y*t

> 

yt := v__y*t;

v__y*t

At time zero, there is only kinetic energy

> 

E__0 := (1/2)*m*(v__0^2+v__y^2);

(1/2)*m*(v__0^2+v__y^2)

At time t, the energy must be the same

> 

eq := E__0 = (1/2)*m*(v__x^2+v__y^2)+m*g*v__y*t

(1/2)*m*(v__0^2+v__y^2) = (1/2)*m*(v__x^2+v__y^2)+m*g*v__y*t

So the x component of the velocity as a function of time is

> 

vt__x := solve(eq, v__x)[1];

(-2*g*t*v__y+v__0^2)^(1/2)

And this will be zero at time t__0, after which the problem doesn't make sense

> 

t__0 := solve(vt__x, t);

(1/2)*v__0^2/(g*v__y)

So the x coordinate is

> 

xt := `assuming`([int(vt__x, t = 0 .. t)], [v__0::positive]);

-(1/3)*((-2*g*t*v__y+v__0^2)^(3/2)-v__0^3)/(g*v__y)

> 

x__0 := eval(xt, t = t__0);

(1/3)*v__0^3/(g*v__y)

And the y coordinate at the end is

> 

y__0 := v__y*t__0;

(1/2)*v__0^2/g

Eliminate t to get y(x)

> 

ans := solve({x = xt, y = yt}, {t, y}, explicit)[1];

{t = -(1/2)*((-3*g*v__y*x+v__0^3)^(2/3)-v__0^2)/(g*v__y), y = -(1/2)*((-3*g*v__y*x+v__0^3)^(2/3)-v__0^2)/g}

This is a shifted form of the semicubical parabola.

> 

simplify(eval((y-y__0)^3/(x-x__0)^2, ans));

-(9/8)*v__y^2/g

> 

NULL

 

Download Niele.mw

Edit: The conclusion is that to get the conventional form, take (x0,y0) at the origin, which means to take v0 =0, as deduced by Rouben. A more straightforward way rather than going backwards is here:

Niele2.mw

coordinateview...

dharr 8270
June 05 2022
3 0

cos(t) goes to zero, meaning 2/cos(t) goes to infinity. So restrict the range with coordinateview. The answer is a vertical line, as you said.
 

> 

plots:-polarplot(2/cos(t), t = 0 .. 2*Pi,coordinateview=[0..6,0..2*Pi],color=red)

> 

``

 

Download polarplot.mw

First 41 42 43 44 45 46 47 Last Page 43 of 82