dharr

Dr. David Harrington

8482 Reputation

22 Badges

21 years, 34 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

Use

fprintf(fd,"%a",u11);
close(fd);

Edit: If you need the "u11:=", then that is

fprintf(fd,"u11:=%a",u11);
close(fd);

I interpret your x1prime as a derivative with respect to time, and assume the partial derivatives don't depend explicitly on time. Then the equations you give aren't really relevant, and something like the following works.

de1:=diff(H(x,y),x)=2*x*y^3;
de2:=diff(H(x,y),y)=3*x^2*y^2;
pdsolve({de1,de2});

{H(x, y) = x^2*y^3+_C1}

Download pdsolve.mw

The inttrans package in my version of Maple (2015) does not work this one out, but it can be worked out step by step using Cauchy's residue theorem.

Edit: redone to use strict upper half plane.

Download Fourier.mw

Here I read the whole file in as a string and then process it (including removing all the spaces). If you really have huge files, then reading line by line as in @mmcdara's answer is superior. But either way, I think the use of sscanf or fscanf helps.

readdata.mw

 

How about

foo:=proc(n::integer);
  if n=0 then
     return FAIL
  else
     return 1,2;
  fi;
end proc:
L:=foo(6);
if L<>FAIL then
  A,B:=L;
end if;

This also works if you return a list and use A,B:=L[]; 

If you really like the call form a,b:=foo(), then you could just return FAIL,FAIL and test A for FAIL.

The nonlinearsimplex method does not use derivatives and so can avoid this problem. But it can't do bounded problems. It finds the maxima at x=0, but these go to infinity, so not sure what you really want here. To find those, just look for where the reciprocal goes to zero, e.g., with fsolve.

plot.mw

restart;
with(GraphTheory):
G:=CycleGraph(4):
SetVertexPositions(G,[[0,0],[1,1],[2,0],[1,-1]]):
G:=RelabelVertices(G,["left","top","right","bottom"]):
DrawGraph(G);
ChromaticPolynomial(G, 5);

260

Download Flag.mw

You can recklessly force this, or play it safe with CancelInverses.

SolveTools:-CancelInverses(arcsin(sin(x)));

x
 

SolveTools:-CancelInverses(arcsin(sin(x)),'safe');

arcsin(sin(x))

 

Not sure I exactly understand how the code in test.mpl is being run. The command currentdir() finds the path of the currently running file (usually a (preexisting) worksheet but perhaps running from a command prompt?), and the subdirectory can then appended.

In my old version 2015, there are no multigraphs or loops, so the code for RandomDigraph produces only simple graphs. Here it is.

DiGraph.mw

I put the PDEs and conditions etc in the correct form, reduced the number of initial conditions to one, (one dt) and the number of boundary conditions to one (one dx). Then it says "Error, (in pdsolve/numeric/match_PDEs_BCs) cannot handle systems with multiple PDE describing the time dependence of the same dependent variable, or having no time dependence". My interpretation is that since the first PDE doesn't have any time derivatives, it can't be handled.

pdsolve_numeric.mw

In my (older) version of Maple simplify(CONV_2) works.
 

CONV_2 := 2*f__2(r, theta, phi)*v__r/(r*sqrt(r^4*sin(phi)^2))+`v__&phi;`*(diff(f__2(r, theta, phi), phi))/sqrt(r^4*sin(phi)^2)+`v__&phi;`*f__2(r, theta, phi)*(diff(1/sqrt(r^4*sin(phi)^2), phi))+cos(phi)*f__2(r, theta, phi)*`v__&phi;`/(sin(phi)*sqrt(r^4*sin(phi)^2))+(diff(f__2(r, theta, phi), theta))*`v__&theta;`/sqrt(r^4*sin(phi)^2)+v__r*(diff(f__2(r, theta, phi), r))/sqrt(r^4*sin(phi)^2)+v__r*f__2(r, theta, phi)*(diff(1/sqrt(r^4*sin(phi)^2), r)):

simplify(CONV_2);

-(cos(phi)^2-1)*(`v__&phi;`*(diff(f__2(r, theta, phi), phi))+(diff(f__2(r, theta, phi), theta))*`v__&theta;`+v__r*(diff(f__2(r, theta, phi), r)))/(sin(phi)^2*(r^4*sin(phi)^2)^(1/2))

``

Download simplify.mw

You had extra {} in eq (5). I also removed the "=0", which is not required and simplifies subsequent manipulations. Then I converted to powers of csc(eta) and lambda only, but it doesn't seem to be in the form you wanted. Note your sqrt(C) is just cot(eta).

Edit: I put back the cot(eta) = sqrt(C) that I had divided out, so it is in the required form if you are OK with B as a polynomial in csc(eta).

restart

with(PDEtools):

declare(u(x, t));

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

(1)

PDE := diff(u(x, t), t)-(diff(u(x, t), `$`(x, 2), t))+3*u(x, t)^2*(diff(u(x, t), x))-2*(diff(u(x, t), x))*(diff(u(x, t), `$`(x, 2)))-u(x, t)*(diff(u(x, t), `$`(x, 3)));

diff(u(x, t), t)-(diff(diff(diff(u(x, t), t), x), x))+3*u(x, t)^2*(diff(u(x, t), x))-2*(diff(u(x, t), x))*(diff(diff(u(x, t), x), x))-u(x, t)*(diff(diff(diff(u(x, t), x), x), x))

(2)

ODE := convert(simplify(eval(PDE, u(x, t) = U(t*lambda__2+x*lambda__1)), {t*lambda__2+x*lambda__1 = eta}), diff);

(-2*(diff(U(eta), eta))*(diff(diff(U(eta), eta), eta))-(diff(diff(diff(U(eta), eta), eta), eta))*U(eta))*lambda__1^3-(diff(diff(diff(U(eta), eta), eta), eta))*lambda__1^2*lambda__2+3*U(eta)^2*(diff(U(eta), eta))*lambda__1+(diff(U(eta), eta))*lambda__2

(3)

ODE1 := simplify((-2*(diff(U(eta), eta))*(diff(diff(U(eta), eta), eta))-(diff(diff(diff(U(eta), eta), eta), eta))*U(eta))*lambda__1^3-(diff(diff(diff(U(eta), eta), eta), eta))*lambda__1^2*lambda__2+3*U(eta)^2*(diff(U(eta), eta))*lambda__1+(diff(U(eta), eta))*lambda__2, 'symbolic')

-2*(diff(U(eta), eta))*(diff(diff(U(eta), eta), eta))*lambda__1^3-(diff(diff(diff(U(eta), eta), eta), eta))*U(eta)*lambda__1^3+3*U(eta)^2*(diff(U(eta), eta))*lambda__1-(diff(diff(diff(U(eta), eta), eta), eta))*lambda__1^2*lambda__2+(diff(U(eta), eta))*lambda__2

(4)

anst := U(eta) = sum(a[i]*csc(eta)^i, i = -2 .. 2);

U(eta) = a[-2]/csc(eta)^2+a[-1]/csc(eta)+a[0]+a[1]*csc(eta)+a[2]*csc(eta)^2

(5)

Eq := normal(subs(anst, ODE1));

-cot(eta)*(-48*csc(eta)^10*cot(eta)^2*lambda__1^3*a[2]^2+15*csc(eta)^11*lambda__1*a[1]*a[2]^2-10*csc(eta)^8*cot(eta)^2*lambda__1^3*a[1]^2+12*csc(eta)^10*lambda__1*a[0]*a[2]^2+12*csc(eta)^10*lambda__1*a[1]^2*a[2]-29*csc(eta)^9*lambda__1^3*a[1]*a[2]+9*csc(eta)^9*lambda__1*a[-1]*a[2]^2-16*csc(eta)^8*lambda__1^3*a[0]*a[2]-16*csc(eta)^8*lambda__1^2*lambda__2*a[2]+6*csc(eta)^8*lambda__1*a[-2]*a[2]^2+6*csc(eta)^8*lambda__1*a[0]^2*a[2]+6*csc(eta)^8*lambda__1*a[0]*a[1]^2-9*csc(eta)^7*lambda__1^3*a[-1]*a[2]-5*csc(eta)^7*lambda__1^3*a[0]*a[1]-5*csc(eta)^7*lambda__1^2*lambda__2*a[1]+3*csc(eta)^7*lambda__1*a[-1]*a[1]^2+3*csc(eta)^7*lambda__1*a[0]^2*a[1]-8*csc(eta)^6*lambda__1^3*a[-2]*a[2]-2*csc(eta)^6*lambda__1^3*a[-1]*a[1]-5*csc(eta)^5*lambda__1^3*a[-2]*a[1]-csc(eta)^5*lambda__1^3*a[-1]*a[0]-csc(eta)^5*lambda__1^2*lambda__2*a[-1]-3*csc(eta)^5*lambda__1*a[-1]^2*a[1]-3*csc(eta)^5*lambda__1*a[-1]*a[0]^2-8*csc(eta)^4*lambda__1^3*a[-2]*a[0]+8*csc(eta)^2*cot(eta)^2*lambda__1^3*a[-2]^2-8*csc(eta)^4*lambda__1^2*lambda__2*a[-2]-6*csc(eta)^4*lambda__1*a[-2]^2*a[2]-6*csc(eta)^4*lambda__1*a[-2]*a[0]^2-6*csc(eta)^4*lambda__1*a[-1]^2*a[0]-17*csc(eta)^3*lambda__1^3*a[-2]*a[-1]-9*csc(eta)^3*lambda__1*a[-2]^2*a[1]-12*csc(eta)^2*lambda__1*a[-2]^2*a[0]-12*csc(eta)^2*lambda__1*a[-2]*a[-1]^2-15*csc(eta)*lambda__1*a[-2]^2*a[-1]-18*csc(eta)^3*lambda__1*a[-2]*a[-1]*a[0]-50*csc(eta)^9*cot(eta)^2*lambda__1^3*a[1]*a[2]-24*csc(eta)^8*cot(eta)^2*lambda__1^3*a[0]*a[2]-24*csc(eta)^8*cot(eta)^2*lambda__1^2*lambda__2*a[2]-12*csc(eta)^7*cot(eta)^2*lambda__1^3*a[-1]*a[2]-6*csc(eta)^7*cot(eta)^2*lambda__1^3*a[0]*a[1]+18*csc(eta)^9*lambda__1*a[0]*a[1]*a[2]-6*csc(eta)^7*cot(eta)^2*lambda__1^2*lambda__2*a[1]-8*csc(eta)^6*cot(eta)^2*lambda__1^3*a[-2]*a[2]-2*csc(eta)^6*cot(eta)^2*lambda__1^3*a[-1]*a[1]+12*csc(eta)^8*lambda__1*a[-1]*a[1]*a[2]-2*csc(eta)^5*cot(eta)^2*lambda__1^3*a[-2]*a[1]+6*csc(eta)^7*lambda__1*a[-2]*a[1]*a[2]+6*csc(eta)^7*lambda__1*a[-1]*a[0]*a[2]+4*csc(eta)^3*cot(eta)^2*lambda__1^3*a[-2]*a[-1]-6*csc(eta)^5*lambda__1*a[-2]*a[-1]*a[2]-6*csc(eta)^5*lambda__1*a[-2]*a[0]*a[1]-12*csc(eta)^4*lambda__1*a[-2]*a[-1]*a[1]-6*lambda__1*a[-2]^3+6*csc(eta)^12*lambda__1*a[2]^3-24*csc(eta)^10*lambda__1^3*a[2]^2+3*csc(eta)^9*lambda__1*a[1]^3-7*csc(eta)^8*lambda__1^3*a[1]^2+2*csc(eta)^8*lambda__2*a[2]+csc(eta)^7*lambda__2*a[1]-3*csc(eta)^4*lambda__1^3*a[-1]^2-csc(eta)^5*lambda__2*a[-1]-3*csc(eta)^3*lambda__1*a[-1]^3-16*csc(eta)^2*lambda__1^3*a[-2]^2-2*csc(eta)^4*lambda__2*a[-2])/csc(eta)^6

(6)

Note that

is(cot(x)^2 = csc(x)^2-1);

true

(7)

So we can convert the cot to csc and leave only powers of csc(eta) and lambda with

cot(eta)*simplify(numer(Eq)/cot(eta), {cot(eta)^2 = csc(eta)^2-1})

cot(eta)*((48*lambda__1^3*a[2]^2-6*lambda__1*a[2]^3)*csc(eta)^12+(50*lambda__1^3*a[1]*a[2]-15*lambda__1*a[1]*a[2]^2)*csc(eta)^11+(24*lambda__1^3*a[0]*a[2]+10*lambda__1^3*a[1]^2-24*lambda__1^3*a[2]^2+24*lambda__1^2*lambda__2*a[2]-12*lambda__1*a[0]*a[2]^2-12*lambda__1*a[1]^2*a[2])*csc(eta)^10+(12*lambda__1^3*a[-1]*a[2]+6*lambda__1^3*a[0]*a[1]-21*lambda__1^3*a[1]*a[2]+6*lambda__1^2*lambda__2*a[1]-9*lambda__1*a[-1]*a[2]^2-18*lambda__1*a[0]*a[1]*a[2]-3*lambda__1*a[1]^3)*csc(eta)^9+(8*lambda__1^3*a[-2]*a[2]+2*lambda__1^3*a[-1]*a[1]-8*lambda__1^3*a[0]*a[2]-3*lambda__1^3*a[1]^2-8*lambda__1^2*lambda__2*a[2]-6*lambda__1*a[-2]*a[2]^2-12*lambda__1*a[-1]*a[1]*a[2]-6*lambda__1*a[0]^2*a[2]-6*lambda__1*a[0]*a[1]^2-2*lambda__2*a[2])*csc(eta)^8+(2*lambda__1^3*a[-2]*a[1]-3*lambda__1^3*a[-1]*a[2]-lambda__1^3*a[0]*a[1]-lambda__1^2*lambda__2*a[1]-6*lambda__1*a[-2]*a[1]*a[2]-6*lambda__1*a[-1]*a[0]*a[2]-3*lambda__1*a[-1]*a[1]^2-3*lambda__1*a[0]^2*a[1]-lambda__2*a[1])*csc(eta)^7+(-4*lambda__1^3*a[-2]*a[-1]+3*lambda__1^3*a[-2]*a[1]+lambda__1^3*a[-1]*a[0]+lambda__1^2*lambda__2*a[-1]+6*lambda__1*a[-2]*a[-1]*a[2]+6*lambda__1*a[-2]*a[0]*a[1]+3*lambda__1*a[-1]^2*a[1]+3*lambda__1*a[-1]*a[0]^2+lambda__2*a[-1])*csc(eta)^5+(-8*lambda__1^3*a[-2]^2+8*lambda__1^3*a[-2]*a[0]+3*lambda__1^3*a[-1]^2+8*lambda__1^2*lambda__2*a[-2]+6*lambda__1*a[-2]^2*a[2]+12*lambda__1*a[-2]*a[-1]*a[1]+6*lambda__1*a[-2]*a[0]^2+6*lambda__1*a[-1]^2*a[0]+2*lambda__2*a[-2])*csc(eta)^4+(21*lambda__1^3*a[-2]*a[-1]+9*lambda__1*a[-2]^2*a[1]+18*lambda__1*a[-2]*a[-1]*a[0]+3*lambda__1*a[-1]^3)*csc(eta)^3+(24*lambda__1^3*a[-2]^2+12*lambda__1*a[-2]^2*a[0]+12*lambda__1*a[-2]*a[-1]^2)*csc(eta)^2+15*csc(eta)*lambda__1*a[-2]^2*a[-1]+6*lambda__1*a[-2]^3)

(8)

Download CM_TW.mw

Many things can be adjusted. (Worksheet not displaying right now).

pts := [1,3,5,9]; n := numelems(pts);
plots:-pointplot(pts, [0$n], symbolsize=20, symbol=solidcircle,
       axis[2]=[tickmarks=0,color=white], color=red, view=[0..10,-0.2..0.2]);

points.mw

 

EIS, my favourite subject! You can just use simplify rather than use a pade approx. To find the elements in terms of the a[i], b[i], simplify,identity usually is the way to go; see ?simplify,identity. (As @acer points out, you don't have to divide through first.) Since there are fewer elements than coefficients, I just used an ad-hoc method that successive removes elements and simplifies the circuit. My sigma and Rw appear different from in your .pdf, but are the same as @acers, so perhaps they are just equivalent forms.

restart

Y := 1/(R__s+1/(s*C__dl+1/(R__ct+1/(sqrt(s)/sigma+1/R__w))))

1/(R__s+1/(s*C__dl+1/(R__ct+1/(s^(1/2)/sigma+1/R__w))))

Let x=sqrt(s)

Y__1 := `assuming`([simplify(eval(Y, s = x^2))], [x::positive]);

(C__dl*R__ct*R__w*x^3+C__dl*R__ct*sigma*x^2+C__dl*R__w*sigma*x^2+R__w*x+sigma)/(C__dl*R__ct*R__s*R__w*x^3+C__dl*R__ct*R__s*sigma*x^2+C__dl*R__s*R__w*sigma*x^2+R__ct*R__w*x+R__s*R__w*x+R__ct*sigma+R__s*sigma+R__w*sigma)

W__1 := add(a[i]*x^i, i = 0 .. 3)/(1+add(b[i]*x^i, i = 1 .. 3));

(x^3*a[3]+x^2*a[2]+x*a[1]+a[0])/(x^3*b[3]+x^2*b[2]+x*b[1]+1)

vars1 := `minus`(indets(Y__1), {x});

{C__dl, R__ct, R__s, R__w, sigma}

{a[0], a[1], a[2], a[3], b[1], b[2], b[3]}

Different number of variables. Easy to get a__i and b__i in terms of the circuit elements

eqs := solve(identity(Y__1 = W__1, x), vars2);

{a[0] = 1/(R__ct+R__s+R__w), a[1] = R__w/(sigma*(R__ct+R__s+R__w)), a[2] = C__dl*(R__ct+R__w)/(R__ct+R__s+R__w), a[3] = C__dl*R__ct*R__w/(sigma*(R__ct+R__s+R__w)), b[1] = R__w*(R__ct+R__s)/(sigma*(R__ct+R__s+R__w)), b[2] = R__s*C__dl*(R__ct+R__w)/(R__ct+R__s+R__w), b[3] = R__s*C__dl*R__ct*R__w/(sigma*(R__ct+R__s+R__w))}

But not the other way around. (Normally this would be the way to find how vars1 relate to vars2 if there was a unique solution.).

solve(identity(Y__1 = W__1, x), vars1)

Let's do it ad-hoc by incrementally extracting the elements.  
Look at the high frequency limit.

limit(Y__1, x = infinity) = limit(W__1, x = infinity);

1/R__s = a[3]/b[3]

R__s = b[3]/a[3]

Subtract off Rs from the impedance, then convert back to admittance

Y__2 := 1/(1/Y-R__s);

s*C__dl+1/(R__ct+1/(s^(1/2)/sigma+1/R__w))

-(x^3*a[3]+x^2*a[2]+x*a[1]+a[0])*a[3]/(x^2*a[2]*b[3]-x^2*a[3]*b[2]+x*a[1]*b[3]-x*a[3]*b[1]+a[0]*b[3]-a[3])

Comparing with the following, we see from the coeff of x^2 in the denominator above that a[2]*b[3]-a[3]*b[2]=0, which gives the second form for Rs

`assuming`([simplify(eval(Y__2, s = x^2))], [positive]);

(C__dl*R__ct*R__w*x^3+C__dl*R__ct*sigma*x^2+C__dl*R__w*sigma*x^2+R__w*x+sigma)/(R__ct*R__w*x+R__ct*sigma+R__w*sigma)

W__2b := simplify(algsubs(a[2]*b[3] = a[3]*b[2], W__2));

-(x^3*a[3]+x^2*a[2]+x*a[1]+a[0])*a[3]/(x*a[1]*b[3]-x*a[3]*b[1]+a[0]*b[3]-a[3])

To find Cdl

Cdleq := limit(Y__2/s, s = infinity) = limit(W__2b/x^2, x = infinity);

C__dl = -a[3]^2/(a[1]*b[3]-a[3]*b[1])

Remove Cdl and convert to impedance

Z__3 := 1/(-C__dl*s+Y__2);

R__ct+1/(s^(1/2)/sigma+1/R__w)

((a[1]*b[3]-a[3]*b[1])*x+a[0]*b[3]-a[3])*(a[1]*b[3]-a[3]*b[1])/(a[3]*((a[0]*a[3]*b[3]-a[1]*a[2]*b[3]+a[2]*a[3]*b[1]-a[3]^2)*x^2+(-a[1]^2*b[3]+a[1]*a[3]*b[1])*x-a[0]*a[1]*b[3]+a[0]*a[3]*b[1]))

Comparing with below, we see that the coefficient of x^2 in the denominator must be zero

`assuming`([simplify(eval(Z__3, s = x^2))], [positive]);

(R__ct*R__w*x+R__ct*sigma+R__w*sigma)/(R__w*x+sigma)

Remove the x^2 term and eliminate a[2]

solve(coeff(denom(W__3), x^2), a[2]);

a[3]*(a[0]*b[3]-a[3])/(a[1]*b[3]-a[3]*b[1])

-(x*a[1]*b[3]-x*a[3]*b[1]+a[0]*b[3]-a[3])/((x*a[1]+a[0])*a[3])

To find Rct

Rcteq := limit(Z__3, s = infinity) = limit(W__3b, x = infinity);

R__ct = -(a[1]*b[3]-a[3]*b[1])/(a[1]*a[3])

Remove Rct and convert to admittance

Y__4 := 1/(Z__3-R__ct);

s^(1/2)/sigma+1/R__w

-a[1]^2*x/(a[0]*b[1]-a[1])-a[0]*a[1]/(a[0]*b[1]-a[1])

Finally, we extract sigma and Rw

sigmaeq := sigma = 1/coeff(W__4, x);

sigma = -(a[0]*b[1]-a[1])/a[1]^2

Rweq := R__w = expand(1/(eval(W__4, x = 0)))

R__w = -b[1]/a[1]+1/a[0]

NULL

Download electrochem.mw

First 38 39 40 41 42 43 44 Last Page 40 of 83