dharr

Dr. David Harrington

6416 Reputation

21 Badges

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

Social Networks and Content at Maplesoft.com

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

Your proposed solution (9) contains the unknown function R(xi), so I'm not sure what you expect here. You can use

simplify(eval(SO, f))

to see the resulting ode in R(xi). I'm not sure if this is the question you are asking.

Since you are using floating point numbers, an alternative is

Optimization:-Minimize(G(x,y), x = 0..1, y = 0..1)

In principle it might find a local rather than global minimum, but a plot shows it is fine in this case.

@Scot Gould 's solution is nice (Vote up), but it is possible to put the Vectors directly in the Matrix if you first make an Array.

restart

vectors1 := [seq(`<,>`(i, i+1, i+2), i = 1 .. 3)]; vectors2 := [seq(`<,>`(i, 2*i, 3*i), i = 1 .. 3)]

[Vector(3, {(1) = 1, (2) = 2, (3) = 3}), Vector(3, {(1) = 2, (2) = 3, (3) = 4}), Vector(3, {(1) = 3, (2) = 4, (3) = 5})]

[Vector[column](%id = 36893491208837358164), Vector[column](%id = 36893491208837358284), Vector[column](%id = 36893491208837358404)]

Arrays accept Vectors as entries

arrowList := Array(1 .. 3, 1 .. 2, [seq([vectors1[n], vectors2[n]], n = 1 .. 3)])

Array(%id = 36893491208837351772)

And we can convert that to a Matrix

arrowList := Matrix(arrowList)

Matrix(%id = 36893491208837333588)

arrowList[3, 2]

Vector[column](%id = 36893491208837358404)

arrowList[3, 2][1]

3

NULL

Download MaplePrimes_Matrix_of_Vectors.mw

The main error here was that psi(x,t) was already assigned when you did pdetest. I don't know the assumptions about the signs of things, but to make progress you will need to make them explicit. I think the presence of the absolute values in pde is going to get you into trouble if you are using complex functions. You had an earlier ode in U(xi) that you derived from the pde. So I would make sure that odetest on that works before attempting the full pde. 

You had kt which I changed to k*t but in your notes it is just k. Isn't all the time dependence in the exponential?

pde-solve.mw

solve, identity will automatically do all the work for you

solve(identity(E1, xi), {k, lambda, w, A[0], A[1], B[1]})

See the attached worksheet for the results.

solve_identity.mw

The pseudo inverse of A seems to work here, but off-hand I can't justify this.

Groebner.mw

 

Data is in the startup code as a hack since you didn't give the datafile. I used fsolve and directly in p__eng1, but NonlinearFit still doesn't return in reasonable time, though you might want to wait longer than I did. As usual in nonlinearfitting, a good initial guess at the solution is probably needed, and reducing the ranges over which you search, if you know them.

test.mw

 

You can remove the solutions you don't want with

remove(has,{COEFFS},[A[0]=0, A[1]=0, B[1]=0]);

If you want to get rid of the RootOfs, add explicit to solve. But if you know the signs of the other parameters add them after the solve, e.g., solve(...) assuming k>0;

eqns.mw

int(F[1](x,y),[x=0..1.,y=0..1.]);

returns 0.2080471649 (This is Maple 2024, which also has a problem with the nested form)

You eq (5) label refers to length(%), not the collected expression, so delete length(%). Change coeff(.., lambda^0) to coeff(..,lambda,0).

But your collect is on a ratio. If you want to just work with the numerator, that can be done, as on the attached.

coeff.mw

The following works, though it is slow.

plots:-inequal({-1 < lambda1, -1 < lambda2, lambda1 < 1, lambda2 < 1}, v = -5 .. 5, z = -5 .. 5)


The colored regions are where both eigenvalues are between -1 and 1. Depending on the ranges you want for v and z, you can refine this, and there are other options for inequal that might increase resolution or efficiency.

2d_implicit_plot_[v_z].mw

 

For FriendshipGraph in SpecialGraphs in GraphTheory GraphTheory[SpecialGraphs][FriendshipGraph] works.

I'm not sure if this is an answer to the specific question, but here is my take on it. Basically, RootOf is very reliable for roots of polynomials, but if there are square roots or cube roots inside it, it is better IMO to reformulate it without these. Here is a way to do that for this case.

dsolve without method

ode:=diff(y(x), x) = (3*x - y(x) + 1)/(3*y(x) - x + 5);
ic:=y(0)=0;
dsolve({ode,ic},implicit);
isol:=eval(lhs(%),y(x)=y);

diff(y(x), x) = (3*x-y(x)+1)/(3*y(x)-x+5)

y(0) = 0

-(2/3)*ln(-(3+y(x)+x)/(x+1))-(1/3)*ln((-1-y(x)+x)/(x+1))-ln(x+1)+(2/3)*ln(3)+I*Pi = 0

-(2/3)*ln(-(3+y+x)/(x+1))-(1/3)*ln((-1-y+x)/(x+1))-ln(x+1)+(2/3)*ln(3)+I*Pi

We need to combine the two arguments of the ln to solve this in the first 2 terms, but the 1/3 powers are a problem

combine(isol,ln);
isol2:=combine(-3*isol,ln);

ln(1/(-(3+y+x)/(x+1))^(2/3))+ln(1/((-1-y+x)/(x+1))^(1/3))-ln((1/3)*(x+1)*3^(1/3))+I*Pi

2*ln(-(3+y+x)/(x+1))+3*ln(x+1)+ln((1/9)*(-1-y+x)/(x+1))-(3*I)*Pi

with_y,no_y:=selectremove(has,isol2,y);

2*ln(-(3+y+x)/(x+1))+ln((1/9)*(-1-y+x)/(x+1)), 3*ln(x+1)-(3*I)*Pi

p1:=simplify(exp(with_y-ln(a)));
p2:=simplify(exp(no_y+ln(a)));

(1/9)*(3+y+x)^2*(-1-y+x)/((x+1)^3*a)

-(x+1)^3*a

So now we have p1*p2=exp(0)

eq:=9*(1-p1*p2);

(-1-y+x)*(3+y+x)^2+9

plots:-implicitplot(eq,x=-10..10,y=-2..9);

solve(.., parametric) suggests that except for the case of x=-1, we have the three solutions of the cubic, the same as we would get just with solve.

solve(eq,y,parametric);

piecewise(x+1 = 0, %SolveTools[Engine]({-y^3-6*y^2-12*y+1}, {y}), x+1 <> 0, [{y = (1/6)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)-(6*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)-(1/3)*x-7/3}, {y = -(1/12)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)+(3*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)-(1/3)*x-7/3+I*sqrt(3)*((1/6)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)+(6*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3))*(1/2)}, {y = -(1/12)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)+(3*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)-(1/3)*x-7/3-I*sqrt(3)*((1/6)*(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3)+(6*(-(4/9)*x^2-(8/9)*x-4/9))/(64*x^3+192*x^2+192*x+1036+36*sqrt(96*x^3+288*x^2+288*x+825))^(1/3))*(1/2)}])

But this seems to miss the real solution for x <~-3 from the plot above - in radical form there are still issues with which square root or which cube root to take

plot([solve(eq,y)],x=-10..10,color=[red, blue, black]);

indexed RootOf is reliable for polynomials

plot([seq(RootOf(eq,y,index=i),i=1..3)],x=-10..10,color=[red, blue, black]);

NULL

 

NULL

Download RootOf_stuff.mw

You may want to add further assumptions.

Download aa.mw

Now that you have corrected the integral with the exp (though one was ecp), the integral is evaluated without any special treatment.

(I had to guess at a value for x.)

restart;

local D; # avoid clash with differential operator

D

Wminus := M*sqrt(Pi)/(2*A*a)*(x*erf(a*x) + (1/(a*sqrt(Pi)))*exp(-a^2*x^2)) + (M/(A*a*sqrt(Pi)))*(Int(exp(-epsilon^2/(4*a^2) - tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2, epsilon = 0..infinity)) + (2*m/Pi)*(Int(exp(-tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2, epsilon = 0..infinity)) + (M/(A*a*sqrt(Pi)))*(Int((1/(2*a^2) + 3*tau*epsilon)*exp(-epsilon^2/(4*a^2) - tau*epsilon^3), epsilon = 0..infinity)) - (2*m/Pi)*tau^(1/3)*GAMMA(2/3) -m*x;

(1/2)*M*Pi^(1/2)*(x*erf(a*x)+exp(-a^2*x^2)/(a*Pi^(1/2)))/(A*a)+M*(Int(exp(-(1/4)*epsilon^2/a^2-tau*epsilon^3)*(cos(epsilon*x)-1)/epsilon^2, epsilon = 0 .. infinity))/(A*a*Pi^(1/2))+2*m*(Int(exp(-tau*epsilon^3)*(cos(epsilon*x)-1)/epsilon^2, epsilon = 0 .. infinity))/Pi+M*(Int(((1/2)/a^2+3*tau*epsilon)*exp(-(1/4)*epsilon^2/a^2-tau*epsilon^3), epsilon = 0 .. infinity))/(A*a*Pi^(1/2))-2*m*tau^(1/3)*GAMMA(2/3)/Pi-m*x

Wplus := M*sqrt(Pi)/(2*A*a)*(x*erf(a*x) + (1/(a*sqrt(Pi)))*exp(-a^2*x^2)) + (M/(A*a*sqrt(Pi)))*evalf(Int(exp(-epsilon^2/(4*a^2) - tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2, epsilon = 0..infinity)) + (2*m/Pi)*evalf(Int(exp(-tau*epsilon^3)*(cos(epsilon*x) - 1)/epsilon^2, epsilon = 0..infinity)) + (M/(A*a*sqrt(Pi)))*evalf(Int((1/(2*a^2) + 3*tau*epsilon)*exp(-epsilon^2/(4*a^2) - tau*epsilon^3), epsilon = 0..infinity)) - (2*m/Pi)*tau^(1/3)*GAMMA(2/3) + m*x;

(1/2)*M*Pi^(1/2)*(x*erf(a*x)+exp(-a^2*x^2)/(a*Pi^(1/2)))/(A*a)+M*(Int(exp(-.2500000000*epsilon^2/a^2-1.*tau*epsilon^3)*(cos(epsilon*x)-1.)/epsilon^2, epsilon = 0. .. Float(infinity)))/(A*a*Pi^(1/2))+2*m*(Int(exp(-1.*tau*epsilon^3)*(cos(epsilon*x)-1.)/epsilon^2, epsilon = 0. .. Float(infinity)))/Pi+M*(Int((.5000000000/a^2+3.*tau*epsilon)*exp(-.2500000000*epsilon^2/a^2-1.*tau*epsilon^3), epsilon = 0. .. Float(infinity)))/(A*a*Pi^(1/2))-2*m*tau^(1/3)*GAMMA(2/3)/Pi+m*x

params:= {x=1e-3, M = 6.3*10^17, t = 7200, a = 10^3, m = 1, D = 1.0*10^(-5), Omega = 1.7*10^23, tau = 1.0*10^(-14)}:
params:= {A = eval(tau/(D*Omega*t),params)} union params;

 

{A = 0.8169934640e-36, D = 0.1000000000e-4, M = 0.6300000000e18, Omega = 0.1700000000e24, a = 1000, m = 1, t = 7200, tau = 0.1000000000e-13, x = 0.1e-2}

evalf(eval(Wminus,params));

0.7711380641e48

evalf(eval(Wplus,params));

0.7711380641e48

NULL

Download integration.mw

2 3 4 5 6 7 8 Last Page 4 of 67