Preben Alsholm

13471 Reputation

22 Badges

20 years, 249 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

Using that the two mixed partials w.r.t. x and y are equal you can write down a condition for solving the system:

restart;
pde1 := diff(f(x,y),x) + A(x) * ( p11* f(x,y) + p12 + p13*x )
= p14 * B(y)/B(x) * diff(f(x,y),x);

pde2 := diff(f(x,y),y) + A(y) * ( p21* f(x,y) + p22 + p23*y )
= p24 * B(x)/B(y) * diff(f(x,y),y);

diff(isolate(pde1,diff(f(x,y),x)),y);
diff(isolate(pde2,diff(f(x,y),y)),x);
factor(rhs(%%-%));
Nm:=numer(%):
#Your first special case:
collect(eval(Nm,{p11=0,p21=0}),diff,factor); #This has to be zero
#Your second special case:
factor(eval(Nm,{p11=0,p21=0,p24=0}));#This has to be zero

In both cases you get requirements on the coefficients for solutions to exist.

In the latter special case the requirement is that

A(x)*B(y)^2*B(x)*(diff(B(y), y))*p14*(p12+p13*x)=0 for all (x,y) in some open non-empty set.

Unless p14 = 0 or p12 = p13 = 0 then the only interesting case would be B = constant.

All parameters must have values, and you need one more initial condition.

I have added D(r)(0)=0.

restart;
G := 1;
M := 1;
pT := 2;

OrbitDEs := diff(r(tau), [tau$2]) = -G*M/r(tau)^2+L^2/r(tau)^3-3*G*M*L^2/r(tau)^4, diff(phi(tau), tau) = L/r(tau)^2, diff(t(tau), tau) = sqrt(r(tau)/(r(tau)-2*G*M)+(diff(r(tau), tau))^2*r(tau)^2/(r(tau)-2*G*M)^2+r(tau)^3*(diff(phi(tau), tau))^2/(r(tau)-2*G*M));

OrbitICs := r(0) = r0, phi(0) = 0, t(0) = 0, D(r)(0)=0;
#Arbitrarily using L = 1, r0=3, and pT= 7.
#You don't have to provide the range argument range = 0..pT.


eqs:=eval({OrbitDEs, OrbitICs},{L=1,r0=3});
OrbitSolution := dsolve(eqs, numeric, {r(tau), phi(tau), t(tau)}, range = 0 .. 7);

plots:-odeplot(OrbitSolution,[tau,phi(tau)],0..3.9);

The error message you reported said to take a look at dsolve,numeric,parameters. That is a good idea, if you intend to change the values of the parameters.

There are two complex solutions to z^2 = a, where a is a given complex number, i.e. two square roots if you will.

Maple's sqrt will find the square root having the least principal argument in absolute value, i.e. the principal square root. (The same is true for other roots).

One square root is sqrt(a) the other is -sqrt(a).

The equation y(lambda,beta) = 0 can be turned into an ode problem

For simplicity:
y:=subs(`βb`=b,y);

eq:=diff(b(lambda),lambda)=-eval(diff(y,lambda)/diff(y,b),b=b(lambda)):

It appears that y is purely imaginary for b between 1 and sqrt(2.1025). For b<1 or b > sqrt(2.1025) and b<sqrt(12.25) y is real.

In the latter two regions the ode can be solved numerically. The initial condition can be found using fsolve with the proper initial point.

I'm using Maple 14.

If I do what you describe, i.e.

s := "T(3)^5*T(5)^35";
t:=parse(s);
u := subs({seq(T(i) = S[i], i = 1 .. 5)}, parse(t));

then I get the error message

Error, invalid input: parse expects its 1st argument, st, to be of type string, but received T(3)^5*T(5)^35

which is reasonable, since t is not a string.

If I replace the last statement by

u := subs({seq(T(i) = S[i], i = 1 .. 5)}, t);

I get the output u := S[3]^5*S[5]^35 as expected.



I'm assuming that you use worksheet interface and 1D-input.
If you intend a statement to take up several lines you should do as you say. No new prompt will appear.

Example:

> for i from 1 to 5 do
   i^2
end do;

Another example:

> plot(    
       sin(x),    
       x=0..Pi    
);

You mixed p and P, and r needs a value.

restart;
with(DEtools):with(plots):

ode:= diff(p(t), t) = r*(1.5-sin(2*t))*p(t);

DEplot(eval(ode,r=1), p(t), t = -3 .. 3, p = 0 .. 1000);

I have never seen that happen. I have seen files where Maple asked the question about how I would like the file opened and I have seen warnings like "there were problems during the loading process, your worksheet may be incomplete", but then there has always been something there when opened in a text editor.

However, the mail system at my university doesn't allow sending mw-files via webmail. By this I mean that if e.g. a student sends me a worksheet by email and I open his mail at home with webmail then the attached file shows the correct size, but there is nothing there, not even in a text editor. When I get to work the next day and open my mail with my mail program , then the file is there and is OK. One solution to this particular problem is to zip the file before sending it.

Can you open the file in a text editor?

If you can somebody might be able to help you. You could upload the file here.

Your last version can be done using unapply, which evaluates its first argument. That helps here (as you noticed).

restart;
with(plots):
Digits:=8;
dPhi:=(rho,z,r,phi)->-1/(1+r^2/2)*r/sqrt(rho^2+z^2+r^2-2*rho*r*cos(phi)):
xmax:=2;xmin:=1;ymax:=2;ymin:=1;gridsize:=5;
f:=(i,j)->Int(Int(diff(dPhi(sqrt(XX[1]^2+XX[2]^2),XX[3],r,phi),XX[i],XX[j]),phi=0..2*Pi),r=0..infinity);
A:=unapply(Matrix(3,f),XX[1],XX[2],XX[3]):
W:=unapply((A(x,y,z)[1,1])^2+(A(x,y,z)[2,2])^2+(A(x,y,z)[3,3])^2+2*(A(x,y,z)[1,2])^2+2*(A(x,y,z)[1,3])^2+2*(A(x,y,z)[2,3])^2,x,y,z):
W2:=unapply(W(x,0,y),x,y):
Vx:=unapply(diff(W2(x,y),x),x,y):
Vy:=unapply(diff(W2(x,y),y),x,y):
V:=(x,y)->[Vx(x,y),Vy(x,y)]:
evalf(V(2,1));
ts:=time():
fieldplot(V(x,y),x=xmin..xmax,y=ymin..ymax,fieldstrength=fixed(0.7),grid=[gridsize,gridsize],colour=grey,axes=boxed,arrows=slim);
ts:=time()-ts:
time=ts;

If you try the following you will see what is going on.

plots:-animate(plot,[[rhs,lhs](A),sigma=-350..5,-130..10],l=3.5..3.6,frames=11);

When you say you are trying to find the solution in the interval [1,10], do you mean sigma or ...?

If you change the sigma interval from -350 .. 5 to-350..-200 you don't see any jump

for K from 3.50 by 0.01  to 3.60 do
    lprint(K, fsolve(subs(l=K, A),sigma=-350..-200));
od;

On the other hand if you do

for K from 3.50 by 0.01  to 3.60 do
>    fsolve(subs(l=K, A),sigma=-200..5);
>    if not has(%,fsolve) then lprint(K,%) end if
> od:

you will get one of the two solutions near -35 (when they exist).

If you know the answer, it can be verified:

A:=arctan(1/(5^(1/2)-1)*(10+2*5^(1/2))^(1/2));
rationalize(tan(A));
convert(tan(2*Pi/5),radical);

Which release of Maple are you using?

I didn't have any problem with Maple 14.01:

restart;
u := 0.4e-1;                           
k := 2.536427;                          
s := 9.466073;
R := 4;
Eq1 := diff(v(t)*sin(theta(t)), t) = k-u*sin(theta(t))*(s*sin(phi(t))+(v(t)*cos(theta(t)))^2/R); Eq2 := diff(v(t)*cos(theta(t)), t) = s*cos(phi(t))-u*cos(theta(t))*(s*sin(phi(t))+(v(t)*cos(theta(t)))^2/R); Eq3 := diff(phi(t), t) = v(t)*cos(theta(t))/R;           
IC := v(0) = 5., theta(0) = .174533, phi(0) = 0.;
vthetafi := dsolve({Eq1, Eq2, Eq3, IC}, {v(t), theta(t), phi(t)}, type = numeric);
vthetafi(1);
plots:-odeplot(vthetafi,[[t,v(t)],[t, theta(t)],[t, phi(t)]],0..10);

It seems easier to use implicitplot:

p:=implicitplot(subs(S0=S01,subs(parvals,A0))=0,l=0..40,R3=-500..0,grid=[50,50]):

L:=op(indets(p,listlist));

Maybe you could use 'events'?

Like this, where there are two levels for S, 40000 and 30000:

piwz2 := piecewise(b(t)=1,(k1-k2)*S(t), k1*S(t));
eq2 := diff(S(t), t) = piwz2;
res:=dsolve({eval(eq2,{k1=1,k2=2}),S(0)=1000,b(0)=0},numeric,events=[[S(t)=40000,b(t)=1],[[S(t)=30000,b(t)=1],b(t)=0]],discrete_variables=[b(t)::boolean]);
plots:-odeplot(res,[t,S(t)],0..10,numpoints=1000);

First 137 138 139 140 141 142 143 Last Page 139 of 158