dharr

Dr. David Harrington

8927 Reputation

22 Badges

21 years, 161 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 replies submitted by dharr

@Andiguys We need to know the sign of the denominator, which means knowing the sign of(8*(Cr*alpha*b-alpha+1)*d+8*alpha*b)*rho0^2+8*(-2*delta+2)*b*d

If delta<1 the second term is positive. For the first term to be positive we need Cr*b>1. If Cr*b is less than 1 then the sign of the first term is harder to decide, there would be more complicated conditions. And if the first term can be negative, we then need to know the relative magnitudes of the first and second terms.

It seems also hard to decide the sign of beta, with the solve f1 introducing new conditions, and the f2 term being too hard for Maple to find conditions.  

So unless you know a lot about the conditions between parameters, I don't think you can get a useful analytical answer. 

@Andiguys The result is correct. I missed the condition t>0, which if you add it reduces the number of options in the output. We have t*(A*t+B) > 0, t>0 and A>0. Since t>0 then A*t+B>0 also, or A*t > -B. Dividing by positive A gives t > -B/A. If B is negative, then -B/A is some positive value and the sole condition is t>-B/A. If B is positive then -B/A is negative; then t>0 and t>-B/A simplifies just to t>0.

Assuming this is the type of analytical result you mean here is the new version. There will be many conditions that need to apply in order to find t. Perhaps you know when some of them apply.

is(dq>=0) gave true on an earlier run for those conditions, but FAIL on rerunning. I think the comment is correct. 

restart

NULL

A := b*(((Cr*alpha*((-g*i2+a)*Cr+2*Crm+2*c-2*Pr)*b+((g*i2-a)*Cr-2*Crm-2*c+2*Pr)*alpha-(-g*i2+a)*Cr)*d+alpha*((-g*i2+a)*Cr+2*Crm+2*c-2*Pr)*b+2*g*i2-2*a)*rho0+(2*((Cr*b-1)*d+b))*(delta+Cn-Pr-1))^2*d/(8*(Cr*b*d+b-d)^2*(((Cr*alpha*b-alpha+1)*d+alpha*b)*rho0^2-2*b*d*(delta-1)))

NULL

B := (Cr*rho0*t*(Cr*alpha*b-alpha-1)*d^2+((alpha*((-g*i2+a)*Cr+2*Crm+2*c+3*t-2*Pr)*Cr*b+((g*i2-a)*Cr-2*Crm-2*c-2*t+2*Pr)*alpha+(g*i2-a)*Cr-2*t)*rho0+(2*(Cr*b-1))*(sigma*t+Cn-Pr+delta-1))*d+(alpha*((-g*i2+a)*Cr+2*Crm+2*c+2*t-2*Pr)*b+2*g*i2-2*a)*rho0+2*b*(sigma*t+Cn-Pr+delta-1))^2*d*b/((8*(((Cr*alpha*b-alpha+1)*rho0^2-2*b*(delta-1))*d+rho0^2*b*alpha))*((Cr*b-1)*d+b)^2)

DATA := [delta = .9, a = 0.1e-1, g = .25, c = 0.5e-1, rho0 = .4, Cn = .4, Crm = .1, i2 = 0.6e-1, alpha = .95, s = 0.1e-1, Pr = .35, upsilon = .95, b = .5, sigma = 0.1e-1, Cr = 0.1e-1, sigma = 0.1e-1, d = .3]

[delta = .9, a = 0.1e-1, g = .25, c = 0.5e-1, rho0 = .4, Cn = .4, Crm = .1, i2 = 0.6e-1, alpha = .95, s = 0.1e-1, Pr = .35, upsilon = .95, b = .5, sigma = 0.1e-1, Cr = 0.1e-1, sigma = 0.1e-1, d = .3]

NULL

We want B-A>=0

q := factor(B-A)

(1/8)*b*d*t*(Cr^2*alpha*b*d^2*rho0+3*Cr*alpha*b*d*rho0-Cr*alpha*d^2*rho0+2*Cr*b*d*sigma-Cr*d^2*rho0+2*alpha*b*rho0-2*alpha*d*rho0+2*b*sigma-2*d*rho0-2*d*sigma)*(Cr^2*alpha*b*d^2*rho0*t-2*Cr^2*alpha*b*d*g*i2*rho0+2*Cr^2*a*alpha*b*d*rho0+4*Cr*Crm*alpha*b*d*rho0-4*Cr*Pr*alpha*b*d*rho0+4*Cr*alpha*b*c*d*rho0+3*Cr*alpha*b*d*rho0*t-2*Cr*alpha*b*g*i2*rho0-Cr*alpha*d^2*rho0*t+2*Cr*alpha*d*g*i2*rho0+2*Cr*a*alpha*b*rho0-2*Cr*a*alpha*d*rho0+2*Cr*b*d*sigma*t-Cr*d^2*rho0*t+2*Cr*d*g*i2*rho0+4*Cn*Cr*b*d-4*Cr*Pr*b*d-2*Cr*a*d*rho0+4*Cr*b*d*delta+4*Crm*alpha*b*rho0-4*Crm*alpha*d*rho0-4*Pr*alpha*b*rho0+4*Pr*alpha*d*rho0+4*alpha*b*c*rho0+2*alpha*b*rho0*t-4*alpha*c*d*rho0-2*alpha*d*rho0*t-4*Cr*b*d+2*b*sigma*t-2*d*rho0*t-2*d*sigma*t+4*g*i2*rho0+4*Cn*b-4*Cn*d-4*Pr*b+4*Pr*d-4*a*rho0+4*b*delta-4*d*delta-4*b+4*d)/((Cr*alpha*b*d*rho0^2+alpha*b*rho0^2-alpha*d*rho0^2-2*b*d*delta+d*rho0^2+2*b*d)*(Cr*b*d+b-d)^2)

Look at the denominator.

dq := simplify(denom(q)); `assuming`([is(dq >= 0)], [positive])

8*((Cr*b-1)*d+b)^2*(((Cr*alpha*rho0^2-2*delta+2)*b-rho0^2*(alpha-1))*d+rho0^2*b*alpha)

FAIL

collect(dq/((Cr*b-1)*d+b)^2, rho0)

(8*(Cr*alpha*b-alpha+1)*d+8*alpha*b)*rho0^2+8*(-2*delta+2)*b*d

So if we knew delta <1 and Cr*b>1 then we would know that dq was nonnegative.

`assuming`([is(dq >= 0)], [positive, delta < 1, Cr*b > 1])

FAIL

So let's assume these conditions apply and the denominator is positive (it can be zero but then the problem doesn't make sense). Then we want the numerator >=0.

nq := collect(numer(q), t, simplify)

d*(Cr*rho0*(Cr*alpha*b-alpha-1)*d^2+((3*Cr*alpha*b-2*alpha-2)*rho0+2*b*Cr*sigma-2*sigma)*d+2*b*(alpha*rho0+sigma))^2*b*t^2+2*(Cr*rho0*(Cr*alpha*b-alpha-1)*d^2+((3*Cr*alpha*b-2*alpha-2)*rho0+2*b*Cr*sigma-2*sigma)*d+2*b*(alpha*rho0+sigma))*d*(((Cr*alpha*((-g*i2+a)*Cr+2*Crm+2*c-2*Pr)*b+((g*i2-a)*Cr-2*Crm-2*c+2*Pr)*alpha-(-g*i2+a)*Cr)*rho0+2*(Cr*b-1)*(delta+Cn-Pr-1))*d+(alpha*((-g*i2+a)*Cr+2*Crm+2*c-2*Pr)*b+2*g*i2-2*a)*rho0+2*b*(delta+Cn-Pr-1))*b*t

It is a quadratic, collect coeffs of t

AA := coeff(nq, t, 2); BB := coeff(nq, t, 1)

d*(Cr*rho0*(Cr*alpha*b-alpha-1)*d^2+((3*Cr*alpha*b-2*alpha-2)*rho0+2*b*Cr*sigma-2*sigma)*d+2*b*(alpha*rho0+sigma))^2*b

2*(Cr*rho0*(Cr*alpha*b-alpha-1)*d^2+((3*Cr*alpha*b-2*alpha-2)*rho0+2*b*Cr*sigma-2*sigma)*d+2*b*(alpha*rho0+sigma))*d*(((Cr*alpha*((-g*i2+a)*Cr+2*Crm+2*c-2*Pr)*b+((g*i2-a)*Cr-2*Crm-2*c+2*Pr)*alpha-(-g*i2+a)*Cr)*rho0+2*(Cr*b-1)*(delta+Cn-Pr-1))*d+(alpha*((-g*i2+a)*Cr+2*Crm+2*c-2*Pr)*b+2*g*i2-2*a)*rho0+2*b*(delta+Cn-Pr-1))*b

`assuming`([is(AA >= 0)], [positive])

true

We don't know the sign of BB

`assuming`([is(BB >= 0)], [positive, delta < 1, Cr*b > 1])

FAIL

We can divide out bd >0 and look at the signs of the remaining factors.

ff := factors(BB/(b*d))

[4, [[Cr^2*alpha*b*d^2*rho0+3*Cr*alpha*b*d*rho0-Cr*alpha*d^2*rho0+2*Cr*b*d*sigma-Cr*d^2*rho0+2*alpha*b*rho0-2*alpha*d*rho0+2*b*sigma-2*d*rho0-2*d*sigma, 1], [-b+d+Cn*b-Cn*d-Pr*b+Pr*d-a*rho0+b*delta-d*delta+(1/2)*Cr^2*a*alpha*b*d*rho0+Cr*Crm*alpha*b*d*rho0-Cr*Pr*alpha*b*d*rho0+Cr*alpha*b*c*d*rho0-(1/2)*Cr*alpha*b*g*i2*rho0+(1/2)*Cr*alpha*d*g*i2*rho0-Cr*b*d+g*i2*rho0+(1/2)*Cr*a*alpha*b*rho0-(1/2)*Cr*a*alpha*d*rho0+(1/2)*Cr*d*g*i2*rho0+Cn*Cr*b*d-Cr*Pr*b*d-(1/2)*Cr*a*d*rho0+Cr*b*d*delta+Crm*alpha*b*rho0-Crm*alpha*d*rho0-Pr*alpha*b*rho0+Pr*alpha*d*rho0+alpha*b*c*rho0-alpha*c*d*rho0-(1/2)*Cr^2*alpha*b*d*g*i2*rho0, 1]]]

f1 := simplify(ff[2][1][1])

Cr*rho0*(Cr*alpha*b-alpha-1)*d^2+((3*Cr*alpha*b-2*alpha-2)*rho0+2*b*Cr*sigma-2*sigma)*d+2*b*(alpha*rho0+sigma)

Is doesn't help (same for <=0)

`assuming`([is(f1 >= 0)], [positive, delta < 1, Cr*b > 1])

FAIL

So we can find some additional conditions - do these apply

`assuming`([solve(f1 > 0, parametric, useassumptions)], [positive, delta < 1, Cr*b > 1])

[[0 < Cr, 0 < d, 1/Cr < b, 0 < rho0, 0 < alpha, alpha < d/(Cr*b*d+b-d), -(1/2)*rho0*(Cr^2*alpha*b*d^2+3*Cr*alpha*b*d-Cr*alpha*d^2-Cr*d^2+2*alpha*b-2*alpha*d-2*d)/(Cr*b*d+b-d) < sigma], [0 < Cr, 0 < d, 1/Cr < b, 0 < rho0, d/(Cr*b*d+b-d) <= alpha, 0 < sigma]]

f2 := ff[2][2][1]

-b+d+Cn*b-Cn*d-Pr*b+Pr*d-a*rho0+b*delta-d*delta+(1/2)*Cr^2*a*alpha*b*d*rho0+Cr*Crm*alpha*b*d*rho0-Cr*Pr*alpha*b*d*rho0+Cr*alpha*b*c*d*rho0-(1/2)*Cr*alpha*b*g*i2*rho0+(1/2)*Cr*alpha*d*g*i2*rho0-Cr*b*d+g*i2*rho0+(1/2)*Cr*a*alpha*b*rho0-(1/2)*Cr*a*alpha*d*rho0+(1/2)*Cr*d*g*i2*rho0+Cn*Cr*b*d-Cr*Pr*b*d-(1/2)*Cr*a*d*rho0+Cr*b*d*delta+Crm*alpha*b*rho0-Crm*alpha*d*rho0-Pr*alpha*b*rho0+Pr*alpha*d*rho0+alpha*b*c*rho0-alpha*c*d*rho0-(1/2)*Cr^2*alpha*b*d*g*i2*rho0

`assuming`([is(f2 >= 0)], [positive, delta < 1, Cr*b > 1])

FAIL

Solving for this factor is time consuming - I didn't wait.

NULL

If we can find all the conditions that make BB positive or negative, then as before (assuming alpha <>0)

solve({alpha > 0, t > 0, t*(alpha*t+beta) > 0}, t, parametric)

piecewise(And(0 < alpha, beta <= 0), [[-beta/alpha < t]], And(0 < alpha, 0 < beta), [[0 < t]], [])

NULL

Download Q_t_Condition.mw

@dharr which leads to a simpler form now given in my answer. Hope this is what you are looking for.

@Andiguys I've corrected the formula at the end and deleted my earlier remark. Since we know that A is positive, then there is no ambiguity; it applies despite the complexity of B or C. So, aside from writing A, B or C in a more simplified form, this seems to be the answer.

@janhardo Thanks. Of course you are right.  And since we knpw the sign of A, there are no other conditions.

@Jean-Michel See my remarks for a possible explanation. I was not answering the question.

For future reference, please upload a worksheet usung the green up-arrow in the Mapleprimes editor.

Your have two parameters sigma1 and sigma2 which don't have values. I'm assuming you want to find these, because you specified two extra boundary conditions. Is that correct, or were these meant to be given values?

Finding parameter values by specifying extra boundary conditions is possible, but much more difficult, often requiring you to give approximate solutions. So you might want to start with specified sigma1 and sigma2, just to make sure you can get some consistent solution.

The error looks a bit like a bug, but when I increase Digits to 20, which bypasses the efficient hardware floating point routines, I find a "matrix is singular" error. This suggests the problem is with the problem setup, not with the solving method. The system is too complicated for me to understand what exactly you are doing. Some boundary conditions look too simple, e.g., setting many values to zero at both boundaries, or some boundary values just proportional to others. Have you set the scale of the problem or does doubling the values of  the c and F quantities lead to another valid solution?

(I am not suggesting you use Digits = 20 other than for debugging purposes.)

You have parameters 100 and 0.01, differing by 4 orders of magnitude, which will make accurate solution difficult. Are the equations non-demensionalized to minimize this issue? Or you could try parameters closer to 1 for an initial solution.

You say it is slow, so does that mean you have a working version that stopped working when you added something. If so, that would be useful information. 

It seems like a tough problem, so perhaps you need to get simpler versions working first (one concentration and potential over [0..1]), if you haven't already done so.

dsolvebvp.mw 

@dharr Ahh, yes, I see how to generalize that:

General.mw

@Harry Garst Your procedure F expects 6 arguments but you gave it 9. Maple ignores those last 3, so taking the first 100 was just taking all combinations of the first 3 and all combinations for the next 3, and can be simplified as:

A4_B5_C4.mw

@Harry Garst For your original question, it seems that the relationship only works for d1 = 3 and taking all the combinations of 2 things out of d2. Trying with random matrices for efficiency  (rather than symbolic ones) for d1=4 and d2=6 there is no value for combinations out of d2 for which the relationship works. I played around a bit, but don't understand why it works in the d1=2; two combination case:

Adjoint353_4.mw

Since you discovered this, I guess you understand better why it works. Since it is not general, it likely doesn't have a name.

Your second example takes the inverse of a block matrix; I don't think this has a special name, but I didn't follow it in detail. Take a look into Schur complements, which may be related.

I understand the 100 possibilities to be 10 x 10, with 10 being the ways to choose 2 things out of 5. Is it 2 because it is d2 - d1? or is it 2 regardless of the size of the matrices.

I don't know of such a formula, but seems bear to  some resemblance to Cauchy-Binet.

Simplified code:

Adjoint353_3.mw

@salim-barzani I don't understand what you mean. I don't understand why the paper doesn't use Eq 2.5 to change Eq 2.8 to an algebraic equation in R and just solve it (for a chosen n). The paper is carelessly prepared and there are some sign issues especially with Eqs 2.2 and 2.5 and the phi equivalent to 2.5 (that @janhardo used, and is not what the paper says). I don't understand about eq 2.2 and the chirp, which just seems to come out of nowhere. Maybe you understand that and can fix the sign issue. Maybe you know how to go from 2.8 to 2.9. This seems like a poor paper to work from to learn this method; I certainly am not spending any more time on it without an explanation of these unresolved issues. Probably just get a different paper that uses the same method and work from that.

@salim-barzani I don't know how they got from 2.8 to 2.9 so I stopped there.

@janhardo I was stuck until I saw your answer, substituting R' instead of phi'. Vote up. Like many similar papers there are lots of typos and poor descriptions.

s1.mw

@salim-barzani The wikipedia page has something about the analogy with trig functions. Its doubly-periodic nature can be shown in a plot:

Doubly periodic nature of the Weierstrass function. For example a hexagonal lattice with unit cell vectors in the complex plane at 60degrees apart. The poles are at the lattice points

See DLMF - 23.5.8 for lattice with omega__1  real, periodicity along the real axis = 2*omega__1. The g-invariants are

restart

g__2 := 0; g__3 := GAMMA(1/3)^18/(4*Pi*`&omega;__1`)^6

0

(64/19683)*Pi^12/(GAMMA(2/3)^18*omega__1^6)

omega_1 = 1.

f := WeierstrassP(z, g__2, eval(g__3, `&omega;__1` = 1))

WeierstrassP(z, 0, (64/19683)*Pi^12/GAMMA(2/3)^18)

p := plots:-complexplot3d(f, z = -5-5*I .. 5+5*I, view = [default, default, 0 .. 4])

Looking down

plots:-display(p, orientation = [-90, 0], style = surfacecontour)

NULL

Download Weierstrass.mw

1 2 3 4 5 6 7 Last Page 1 of 96