## dharr

Dr. David Harrington

## 6416 Reputation

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.

## unknown function...

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.

## Minimize...

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.

## Array, then Matrix...

@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.

 >
 >

Arrays accept Vectors as entries

 >

And we can convert that to a Matrix

 >

 >

 >

 >

## psi assigned...

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...

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

## pseudo-inverse...

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

Groebner.mw

## progress...

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

## remove...

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

## double integral - alternate form...

`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)

## couple of things...

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

## plots:-inequal...

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.

## RootOf reliability...

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);

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);

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

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

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

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

 > 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);

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]);

 >

## assuming...

You may want to add further assumptions.

## Straightforward...

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

 > 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;

 > 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;

 > 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;
 >

 > evalf(eval(Wminus,params));

 > evalf(eval(Wplus,params));

 >