Advanced Search

dharr

Dr. David Harrington

8205 Reputation

22 Badges

20 years, 337 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

treat rows and cols separately...

dharr 8205
August 06 2019
1 1

Like this?


 

> 

restart
> 

rowX, colX := [2, 4], [11, 11, 12, 14, 14, 15]
> 

rowY, colY := [1, 3], [2, 5, 8, 2, 5, 8]
> 

rowlabels, collabels := [A, B], [a, b, c, d, e, f]
> 

m, n := numelems(rowX), numelems(colX)

2, 6

> 

T := Matrix(m, n, proc (i, j) options operator, arrow; sqrt((rowX[i]-colX[j])^2+(rowY[i]-colY[j])^2) end proc)
> 

Distance := DataFrame(T, 'rows' = rowlabels, 'columns' = collabels)

Distance := Matrix(3, 7, {(1, 1) = ``, (1, 2) = a, (1, 3) = b, (1, 4) = c, (1, 5) = d, (1, 6) = e, (1, 7) = f, (2, 1) = A, (2, 2) = 82^(1/2), (2, 3) = 97^(1/2), (2, 4) = 149^(1/2), (2, 5) = 145^(1/2), (2, 6) = 4*10^(1/2), (2, 7) = 218^(1/2), (3, 1) = B, (3, 2) = 5*2^(1/2), (3, 3) = 53^(1/2), (3, 4) = 89^(1/2), (3, 5) = 101^(1/2), (3, 6) = 2*26^(1/2), (3, 7) = 146^(1/2)})

> 

 


 

Download Datafrome.mw

GraphTheory way...

dharr 8205
August 02 2019
0 1

Just for another variation. At least in Maple 2017, I can't export the DataFrame to Excel, even though it can import into a Dataframe.
 

> 

restart;
> 

with(GraphTheory):
> 

X := [3, 6, 4, 13]:   Y := [8, 6, 5, 7]: lbls:=[A,B,C,D];
Pts:=zip((x,y)->[x,y],X,Y);
Pts2:=map(x->x[1]+I*x[2],Pts);
n:=numelems(Pts);

[A, B, C, D]

[[3, 8], [6, 6], [4, 5], [13, 7]]

[3+8*I, 6+6*I, 4+5*I, 13+7*I]

4

> 

dists:=Matrix(n,n,(i,j)->abs(Pts2[i]-Pts2[j]),shape=symmetric):
df:=DataFrame(dists,'rows'=lbls,'columns'=lbls);

df := Matrix(5, 5, {(1, 1) = ``, (1, 2) = A, (1, 3) = B, (1, 4) = C, (1, 5) = D, (2, 1) = A, (2, 2) = 0, (2, 3) = 13^(1/2), (2, 4) = 10^(1/2), (2, 5) = 101^(1/2), (3, 1) = B, (3, 2) = 13^(1/2), (3, 3) = 0, (3, 4) = 5^(1/2), (3, 5) = 5*2^(1/2), (4, 1) = C, (4, 2) = 10^(1/2), (4, 3) = 5^(1/2), (4, 4) = 0, (4, 5) = 85^(1/2), (5, 1) = D, (5, 2) = 101^(1/2), (5, 3) = 5*2^(1/2), (5, 4) = 85^(1/2), (5, 5) = 0})

> 

G := Graph(dists,lbls):
SetVertexPositions(G,Pts):
DrawGraph(G);

> 

 


 

Download Points.mw

An example...

dharr 8205
August 01 2019
0 0

An example of solving a PDE with Laplace transforms (with both initial and boundary conditions) is given in the Application Center:

https://www.maplesoft.com/applications/view.aspx?SID=96958 

improved assumptions...

dharr 8205
July 25 2019
0 1

Try simplify(my_sol) assuming cos(2*x)>=0; Now it really is simpler. But is that valid for your case?

Maple doesn't want to make assumptions that aren't generally true. So if you want to know really what the conditions are, you need to think about why csgn or signum appeared and then make the appropriate assumptions, if they are valid for your problem. They often arise from sqrts of negative arguments, as here. If you want to be more reckless (I usually do this in these cases), just use simplify( , symbolic) and bear in mind that the answer has some unwritten assumptions.

separation of variables solution...

dharr 8205
July 25 2019
1 1

In Maple 2017, the steps through to the standard separation of variable solution can be done step by step. Solution is complex (imaginary)


 

> 

restart;
> 

pde := diff(c(x, t), x, x) - h*diff(c(x, t), x) - diff(c(x, t), t);
iv := c(0, t) = 0, c(a, t) = 0, c(x, 0) = c0;

diff(diff(c(x, t), x), x)-h*(diff(c(x, t), x))-(diff(c(x, t), t))

c(0, t) = 0, c(a, t) = 0, c(x, 0) = c0

Maple suggests a separation of variables.

> 

gensol:=pdsolve(pde);
de1:=op(2,gensol)[][1];

PDESolStruc(c(x, t) = _F1(x)*_F2(t), [{diff(_F2(t), t) = _F2(t)*_c[1], diff(diff(_F1(x), x), x) = _c[1]*_F1(x)+(diff(_F1(x), x))*h}])

diff(_F2(t), t) = _F2(t)*_c[1]

solution of the first de gives an exponential. Presumably t is time and so _c[1] will be negative - call it -lambda

> 

dsolve(de1);

_F2(t) = _C1*exp(_c[1]*t)

and seek solution f(x)*exp(-lambda*t). Just the same as _F1(t) above

> 

de2:=simplify(eval(pde,c(x,t)=f(x)*exp(-lambda*t))/exp(-lambda*t));

-h*(diff(f(x), x))+f(x)*lambda+diff(diff(f(x), x), x)

Solve this de with bc at x=0 c(0,t)=0 => f(0)=0. Divide by the scale factor _C2, which we will put in later (otherwise solving for zero will give the trivial solution _C2=0).

> 

ans:=simplify(rhs(dsolve({de2,f(0)=0}))/_C2);

-exp((1/2)*(h+(h^2-4*lambda)^(1/2))*x)+exp(-(1/2)*(-h+(h^2-4*lambda)^(1/2))*x)

Now we need to find the eigenvalues using the other boundary condition - turns out there is an exact solution. _Z1 means any integer - rename it n, and substitute it back into above.

> 

expand(solve(eval(ans,x=a),lambda,allsolutions));
lambda=subs(_Z1=n,%);
cc:=eval(ans*exp(-lambda*t),%);

(1/4)*h^2+Pi^2*_Z1^2/a^2

lambda = (1/4)*h^2+Pi^2*n^2/a^2

(-exp((1/2)*(h+(-4*Pi^2*n^2/a^2)^(1/2))*x)+exp(-(1/2)*(-h+(-4*Pi^2*n^2/a^2)^(1/2))*x))*exp(-((1/4)*h^2+Pi^2*n^2/a^2)*t)

Solution is complex. Since n=0 gives solution everywhere zero, eigenvalues start with n=1

> 

ccc:=simplify(evalc(cc)) assuming n::posint,a::positive;
eval(ccc,n=0);

-(2*I)*exp(-(1/4)*(4*Pi^2*n^2*t+a^2*h^2*t-2*a^2*h*x)/a^2)*sin(n*Pi*x/a)

0

Check that the solution works for all positive n (except for the initial condition)

> 

map(simplify,pdetest(c(x,t)=ccc,[pde,iv[1..2]])) assuming n::posint,a::positive;

[0, 0, 0]

General solution is a sum of these, where scale factor b is determined by the initial condition (and isn't easily found symbolically)

> 

sol:=c(x,t)=b*Sum(ccc,n=1..infinity);
value(%);

c(x, t) = b*(Sum(-(2*I)*exp(-(1/4)*(4*Pi^2*n^2*t+a^2*h^2*t-2*a^2*h*x)/a^2)*sin(n*Pi*x/a), n = 1 .. infinity))

c(x, t) = b*(sum(-(2*I)*exp(-(1/4)*(4*Pi^2*n^2*t+a^2*h^2*t-2*a^2*h*x)/a^2)*sin(n*Pi*x/a), n = 1 .. infinity))

> 

 


 

Download pdsolve.mw

use Heaviside...

dharr 8205
July 24 2019
0 8

[Edit: incorrect as Carl is pointing out.I put Heaviside(x) instead of x inside the exponential, since exp(-x) goes to infinity at x=-infinity.] Redid using laplace transform and then replaced s=I*omega. Mostly gives what you expect.

Laplace.mw

PDF, mean and stdev...

dharr 8205
July 19 2019
0 4

Like this?
 

> 

X := Statistics:-RandomVariable(Normal(0, 1)) ;

_R

> 

Statistics:-PDF(X,t);
Statistics:-Mean(X);
Statistics:-StandardDeviation(X);

(1/2)*2^(1/2)*exp(-(1/2)*t^2)/Pi^(1/2)

0

1

> 

 


 

Download RandomVariable.mw

DFT...

dharr 8205
July 11 2019
1 2

Fourier series usually means the infinite series, yours looks more like a discrete fourier transform. Search DFT - there is a version in the SignalProcessing package called DFT (or FFT for efficiency) and another in the DiscreteTransforms package called FourierTransform. The formulations are not exactly as you have, but run from 0 to N-1, and have the exponential negative for the forward transform and positive for the inverse transform.

zero...

dharr 8205
July 10 2019
3 2

If I leave tau undefined and change "evalf" to "simplify" so there are no floats, the integral evaluates to zero. Is there a reason to thnk this value is wrong?

070919_ThetaIntegral.mw

f1,1 and f11...

dharr 8205
July 10 2019
1 0

The f_1,1 in the first line is diiferent from the f_11 later on. If you make them all f_11 it works.

2-D notation problem?...

dharr 8205
July 09 2019
1 0

Since you didn't upload your worksheet, it is hard to diagnose, but it works here in 1-D:
 

> 

restart;
> 

with(PDEtools):
> 

declare(u(x,t));
U:=diff_table(u(x,t));

u(x, t)*`will now be displayed as`*u

table( [( ) = u(x, t) ] )

> 

PDE := U[x,x]-U[t];

diff(diff(u(x, t), x), x)-(diff(u(x, t), t))

> 

Infinitesimals(PDE);

[_xi[x](x, t, u) = 0, _xi[t](x, t, u) = 1, _eta[u](x, t, u) = 0], [_xi[x](x, t, u) = 1, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = 0], [_xi[x](x, t, u) = 0, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = u], [_xi[x](x, t, u) = (1/2)*x, _xi[t](x, t, u) = t, _eta[u](x, t, u) = 0], [_xi[x](x, t, u) = -2*t, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = x*u], [_xi[x](x, t, u) = 0, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = exp(_c[1]*t)*exp(_c[1]^(1/2)*x)], [_xi[x](x, t, u) = 0, _xi[t](x, t, u) = 0, _eta[u](x, t, u) = exp(_c[1]*t)/exp(_c[1]^(1/2)*x)], [_xi[x](x, t, u) = (1/2)*x*t, _xi[t](x, t, u) = (1/2)*t^2, _eta[u](x, t, u) = -(1/8)*x^2*u-(1/4)*u*t]

> 

 


 

Download Infinitesimals.mw

identify...

dharr 8205
June 17 2019
1 0

identify(1.77245); gives sqrt(Pi)

asympt of RootOf...

dharr 8205
June 16 2019
2 5
> 

restart;
> 

eq:=((2*n-1)*x + 2*n + 1)*x^n = 1 - x;

((2*n-1)*x+2*n+1)*x^n = 1-x

> 

solve(eq,x);
asympt(%,n);
ass:=convert(%,polynom);

RootOf(2*_Z^n*n*_Z+2*_Z^n*n-_Z^n*_Z+_Z^n+_Z-1)

1-LambertW(4*n^2)/n+(1/2)*LambertW(4*n^2)^2/n^2-(1/12)*LambertW(4*n^2)^2*(2*LambertW(4*n^2)^2+3*LambertW(4*n^2)+3)/((LambertW(4*n^2)+1)*n^3)+(1/24)*LambertW(4*n^2)^3*(LambertW(4*n^2)^2+3*LambertW(4*n^2)+6)/((LambertW(4*n^2)+1)*n^4)-(1/1440)*LambertW(4*n^2)^3*(12*LambertW(4*n^2)^5+89*LambertW(4*n^2)^4+312*LambertW(4*n^2)^3+435*LambertW(4*n^2)^2+270*LambertW(4*n^2)+90)/((LambertW(4*n^2)+1)^3*n^5)+O(1/n^6)

1-LambertW(4*n^2)/n+(1/2)*LambertW(4*n^2)^2/n^2-(1/12)*LambertW(4*n^2)^2*(2*LambertW(4*n^2)^2+3*LambertW(4*n^2)+3)/((LambertW(4*n^2)+1)*n^3)+(1/24)*LambertW(4*n^2)^3*(LambertW(4*n^2)^2+3*LambertW(4*n^2)+6)/((LambertW(4*n^2)+1)*n^4)-(1/1440)*LambertW(4*n^2)^3*(12*LambertW(4*n^2)^5+89*LambertW(4*n^2)^4+312*LambertW(4*n^2)^3+435*LambertW(4*n^2)^2+270*LambertW(4*n^2)+90)/((LambertW(4*n^2)+1)^3*n^5)

> 

n:=1000;
fsolve(eq,x=0..1);
evalf(ass);

1000

.9874167130

.9874167131

> 

 


 

Download limit.mw

ad hoc solution...

dharr 8205
June 15 2019
1 1

If I understand correctly,you know B,and you only know the eigenvalues relative to the most negative one. If so, then you can work out V3 and V6 from the symbolic solutions for the eigenvalues - I'm also assuming you know which are degenerate. Your matrix has a lot of structure, so it might be possible to generalize this.

Download HinderedRotor2.mw

missing multiplication and missing term...

dharr 8205
June 08 2019
2 6

You do not have a space between n and (n+1). Insert a space or multiplication sign between them and you will get a different solution, presumable the one you want. However, (see ?LegendreP), your differential equation needs another term to qualify - adding that leads to a solution in terms of Legendre functions.


 

NULL

NULL

> 

restart
> 

ode := (-x^2+1)*(diff(y(x), x, x))+n*(n+1)*y(x) = 0

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

> 

sol[1] := rhs(dsolve(ode))

_C1*(-x^2+1)*hypergeom([1+(1/2)*n, 1/2-(1/2)*n], [1/2], x^2)+_C2*(-x^3+x)*hypergeom([1-(1/2)*n, 3/2+(1/2)*n], [3/2], x^2)

> 

`assuming`([convert(sol[1], StandardFunctions)], [n::posint])

_C1*(-x^2+1)*hypergeom([1+(1/2)*n, 1/2-(1/2)*n], [1/2], x^2)+_C2*(-x^3+x)*hypergeom([1-(1/2)*n, 3/2+(1/2)*n], [3/2], x^2)

> 

ode2 := (-x^2+1)*(diff(y(x), x, x))-2*x*(diff(y(x), x))+n*(n+1)*y(x) = 0

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

> 

sol[2] := rhs(dsolve(ode2))

_C1*LegendreP(n, x)+_C2*LegendreQ(n, x)

NULL


 

Download convert-Legendre.mw

First 61 62 63 64 65 66 67 Last Page 63 of 81