dharr

Dr. David Harrington

8502 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 replies submitted by dharr

@salim-barzani If I understand you correctly, you are asking where the transformation equations for xp, yp come from.

restart

with(PDEtools)

undeclare(prime, quiet); declare(u(x, y, t), quiet); declare(f(x, y, t), quiet)

pde := diff(u(x, y, t), t)-(diff(diff(u(x, y, t), `$`(x, 4))+5*u(x, y, t)*(diff(u(x, y, t), `$`(x, 2)))+(5/3)*u(x, y, t)^3+5*(diff(u(x, y, t), x, y)), x))-5*u(x, y, t)*(diff(u(x, y, t), y))+5*(int(diff(u(x, y, t), `$`(y, 2)), x))-5*(diff(u(x, y, t), x))*(int(diff(u(x, y, t), y), x))

pde_linear, pde_nonlinear := selectremove(proc (term) options operator, arrow; not has((eval(term, u(x, y, t) = a*u(x, y, t)))/a, a) end proc, expand(pde))

thetai := t*w[i]+y*p[i]+x

eqw := w[i] = -5*p[i]^2

Bij := proc (i, j) options operator, arrow; (-6*p[i]-6*p[j])/(p[i]-p[j])^2 end proc

NULL

theta1 := normal(eval(eval(thetai, eqw), i = 1)); theta2 := normal(eval(eval(thetai, eqw), i = 2))

eqf := f(x, y, t) = theta1*theta2+Bij(1, 2)

eqfcomplex := collect(evalc(eval(eval(eqf, p[2] = conjugate(p[1])), p[1] = a+I*b)), t)

eq17 := u(x, y, t) = 6*(diff(diff(f(x, y, t), x), x))/f(x, y, t)-6*(diff(f(x, y, t), x))^2/f(x, y, t)^2; equ := simplify(eval(eq17, eqfcomplex))

u(x, y, t) = -(12/25)*(b^6*t^2+(-6*a^2*t^2+(6/5)*a*t*y+(2/5)*t*x-(1/25)*y^2)*b^4+(t*a^2-(1/5)*y*a-(1/5)*x)^2*b^2-(3/25)*a)*b^2/(b^6*t^2+(2*a^2*t^2-(2/5)*a*t*y+(2/5)*t*x+(1/25)*y^2)*b^4+(t*a^2-(1/5)*y*a-(1/5)*x)^2*b^2+(3/25)*a)^2

So we want to find a substitution that removes the time dependence from u. One way is to find the maximum and see how it moves. Here, the first solution gives what we want.

ans := solve({diff(rhs(equ), x), diff(rhs(equ), y)}, {x, y})

{x = -5*a^2*t-5*b^2*t, y = 10*a*t}, {x = RootOf(25*a^4*b^2*t^2+50*a^2*b^4*t^2+25*b^6*t^2-9*a+(10*a^2*b*t+10*b^3*t)*_Z+_Z^2)/b, y = 10*a*t}

NULL

NULL

Download Y-pde-line2.mw

@Alfred_F In general, Maple is not going to give information that will verify an output (there are other computer systems for proofs in which every step is justified). Setting infolevel may give information, as @nm showed for the integral in the other thread. Here infolevel in the one-line limit that @Scot Gould showed gives no information. As he showed, you can also look into the actual code used, though this is a painful and often fruitless endeavor.

But for limits one can imagine that asymptotic series is likely to be a useful method. So one could do it partly by hand as below. Of course, then can then ask how Maple got the series, but you probably know the mechanics of this.

restart

Not sure why we have to do this twice.

`assuming`([asympt(x^(1-epsilon)*(int(sin(t^2), t = x .. x+1)), x)], [epsilon > 0])

Error, (in assuming) when calling 'asympt'. Received: 'unable to compute series'

q := `assuming`([asympt(x^(1-epsilon)*(int(sin(t^2), t = x .. x+1)), x)], [epsilon > 0])

(((1/2)*cos(x^2)-(1/2)*cos((x+1)^2))/x+(1/2)*cos((x+1)^2)/x^2+((1/4)*sin(x^2)-(1/4)*sin((x+1)^2)-(1/2)*cos((x+1)^2))/x^3+((3/4)*sin((x+1)^2)+(1/2)*cos((x+1)^2))/x^4+(-(3/8)*cos(x^2)-(1/8)*cos((x+1)^2)-(3/2)*sin((x+1)^2))/x^5+O(1/x^6))*(1/x)^(-1+epsilon)

Throw away terms with larger powers in the denominator, and we see bounded numerator and denominator of x^epsilon -> infinity.

q2 := op([1, 1], q)*op(2, q); simplify(%)

((1/2)*cos(x^2)-(1/2)*cos((x+1)^2))*(1/x)^(-1+epsilon)/x

(1/2)*(cos(x^2)-cos((x+1)^2))*(1/x)^epsilon

NULL

Download limit.mw

@mmcdara Yes, the linear code is fragile, which I solved by expand, so it is still rather fragile; your solution is better. But then the solve part gives the two solutions [0, 2], which doesn't agree with yours, so I'm not sure what is going on.

@salim-barzani The eq15 in the paper seems to work, together with the indefinite integral.

Edit: nonlinear problem pointed out by @mmcdara fixed.

pdeint.mw

I don't think the integral from -infinity to x is any different from your indefinite integral earlier - I would stick to the indefinite integral. Are you sure the ansatz in eq 17 applies to this case?

@acer I missed that because checking the structured types page, i see identical(expr) and not identical(exprseq). After you pointed it out I see it further down in the notes (and vaguely recall using it earlier). I have a expectation that the top of the page will have all calling sequences, but this is not always the case. Makes me miss BNF in help pages (no, not really).

@Christopher2222 Interesting. I didn't think 10 was a type.

@salim-barzani As I said before "I do not yet know how to arrive at the limiting form for fN." That means I tried it but couldn't figure it out. (I don't know how to do it by hand.)

As for when to use simplify, my rule would be:  when you think your answer should be simpler, try simplify.

@mary120 I cannot read the equations and conditions properly. But you said "For some values of parameteres, the output curves were very similar to what I expected". I took that to mean you have some suitable ouput. I would like to see that, because it would clarify the questions I asked and your reply does not help.

Did you try the code provided here? (Hint: you may have to simplify the result.)

@mary120 @mmcdara has suggested you were really wanting a different differential equation than the one you specified in your worksheet, and since you have changed his answer to best I'm assuming that is correct? But now you are saying you can get the correct answers with my worksheet, which uses one branch of the original ode. So is that the right ode and the right branch? Which was the correct initial condition? I suggest you upload a corrected worksheet in which your parameters that do give close to the required solution are given, so we can see the answers to these questions.

The answer to your immediate question is that to get V(phi) values down to phi=0, you need to solve the de all the way out to large enough x values that phi is close to zero, i.e., choose a wider range in

 

x_vals := [seq(i*0.01, i = 0 .. 100)];

From your graph of phi vs x at initial condition phi(0) = 0.022, it looks as though going to x = 250 is about right.

 

@mmcdara You are pobably right, but the OP titled the question with "first oder ode", which I took at face value.

@salim-barzani OK, I understand. Here's how to automate the substitutions. But it is not a solution. Each of your pdes needs some individual attention, to find independent equations. For most cases, the expx in f are the same exps that survive after substituting into the pde, but if thay get messed up, then one needs to work harder, as here. I don't know the general solution to this issue. It may be that simplification with side relations can work; I'll give this some more thought.

Dr.D-pde.mw

I guess my alien powers have deserted me!

@salim-barzani If you only want the variables, and not the more complicated expression, use indets(something, name).

I think this is what you are asking, but if it is not, please be more specific about what you want.

@mmcdara I'm not sure of the point you are making. If it is just that you have to be careful of branch cut issues, I agree. In this case the solve strategy I used to answer the OP works - if both solutions had been real there would have been less of a problem.

(Strangely, the worksheet is not displaying, in a way I haven't seen before.)

Download DE_solve.mw

First 12 13 14 15 16 17 18 Last Page 14 of 89