Preben Alsholm

11719 Reputation

22 Badges

15 years, 254 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

simplify has made use of the form simplify(expr, assume= ...)  before assuming was ever introduced.

They don't work quite the same:
 

restart;
f:=x->a*x;
simplify(sqrt(f(a)),assume=[a<0]);
simplify(sqrt(f(a))) assuming a<0;

The assume = [a<0] version works because simplify evaluates its arguments before it starts its general business as is the normal behavior in Maple.
The assuming version doesn't work here because it really works (as Carl points out) like this:
 

`assuming`([simplify(sqrt(f(a)))],[a<0]);

`assuming`  doesn't evaluate its first argument until the assumption a<0 is made on the visible a, i.e. the other a resulting from evaluating f(a) is not seen at that time, so internally simplify will then be looking at sqrt( a*a~), which it cannot do anything with since a and a~ are different at that time. At exit a~ is replaced by a, and you just have sqrt(a^2).

To use dsolve/numeric/bvp for your new problem you can use the symmetry.
 

restart;
t := -1:
S := 1:
R := 9:

EQ:={(diff(F(x), x $ 4)) + R*( diff(F(x),x)*diff(F(x),x$2)- F(x)*diff(F(x),x$3) ) + R*t*( F(x)*diff(F(x),x$5) - diff(F(x),x)*diff(F(x),x$4) ) - R*S^2*diff(F(x),x$2) =0};


IC:={D(F)(-1)=0, D(F)(1)=0,F(-1)=-1,F(1)=1,D(D(F))(1)=0};

(*Since it follows from EQ and your proposed boundary conditions that there is (could be) an odd solution F, i.e. one that satisfies
F(x) = -F(-x), you can restrict yourself to the interval -1..0. For an odd function F we have F(0) = 0, (D@@2)(F)(0)=0,(D@@4)(F)(0)=0, so we take: *)
 
bcs1:={D(F)(-1)=0,F(0)=0,F(-1)=-1,(D@@2)(F)(0)=0,(D@@4)(F)(0)=0};
sol1:=dsolve(EQ union bcs1,numeric,method=bvp[midrich]);
plots:-odeplot(sol1,[x,F(x)]);
A1:=Array([-1,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0]);
solA1:=dsolve(EQ union bcs1, numeric, method=bvp[midrich],output=A1);
M1:=solA1[2,1];
solA1[1,1]
M2:=Matrix(10,6,(i,j)->`if`(j=1,-M1[11-i,1],(-1)^(j+1)*M1[11-i,j]));
M:=< M1,M2>;
plot(M[..,1..2]);

You can also use odeplot, but then you need to make a structure like solA1:
 

sol:=Matrix(2,1);
sol[1,1]:=solA1[1,1];
sol[2,1]:=M;
plots:-odeplot(sol,[x,F(x)],-1..1);

 

You can do as follows. The names chosen by me can obviously be changed to your liking:
 

restart;
ode1:=diff(u(t), t) = a*u(t) - alpha__1 *u(t)*v(t)/(1 + k*u(t));
ode2:=diff(v(t), t) = -b*v(t) + alpha__1*u(t)*v(t)/(k*u(t) + 1) - alpha__2*v(t)*w(t);
ode3:=diff(w(t), t) = -c*(w(t) - wstar) + alpha__2*v(t)*w(t);
sys1:=PDEtools:-dchange({u(t)=U*p(tau),v(t)=V*q(tau),w(t)=W*r(tau),t=T*tau},{ode1,ode2,ode3},[p,q,r,tau]);
sys2:=solve(sys1,diff~({p,q,r}(tau),tau));
sys3:=expand(sys2);
# Now we pick T, U, and V.
# Looking at the first term on the right of the first ode in sys3 we see from the denominator that U must be chosen so k*U is dimensionless. 
# An obvious choice is U=1/k.
# After that T is chosen as 1/a. Then that first term is handled and is dimensionless.
# Moving to the next term in the same ode we choose W=a/alpha__2:
sys4:=eval(sys3,{U=1/k,T=1/a,V=a/alpha__1,W=a/alpha__2});
# Next step: Find the dimensionless groups and name them:
G:=[b/a,alpha__1/(a*k),c/a,alpha__2*c*wstar/a^2,alpha__2/alpha__1]=~[beta,delta,eta,zeta,a21];
SBS:=solve(G,{a,b,k,c,alpha__1});
## The dimensionless system with the dimensionless parameters beta,delta,eta,zeta,a21:
sys:=subs(SBS,sys4);

The system found is:
sys := {diff(p(tau), tau) = p(tau)^2/(1 + p(tau)) - p(tau)*q(tau)/(1 + p(tau)) + p(tau)/(1 + p(tau)), diff(q(tau), tau) = -q(tau)*r(tau)*p(tau)/(1 + p(tau)) - q(tau)*p(tau)*beta/(1 + p(tau)) - q(tau)*r(tau)/(1 + p(tau)) + q(tau)*p(tau)*delta/(1 + p(tau)) - q(tau)*beta/(1 + p(tau)), diff(r(tau), tau) = a21*q(tau)*r(tau) - eta*r(tau) + zeta}

Separate the text reion from the input region with F3.

Put the cursor in front of de:= ...  and press F3.
After that you could select the whole text line and convert to plain text.
In fact this is all you need to do.

Increase maxmesh. In this case maxmesh = 256 is sufficient.
 

sol:= dsolve(EQ union IC,numeric,maxmesh=256,output=Array([-1,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]));

 

The sum to infinity can be made finite for plotting purposes:
 

restart;

pde := diff(T(x, y), x $ 2) + diff(T(x, y), y $ 2) = 0;
bc := T(0, y) = T1, T(a, y) = T2, T(x, 0) = T2, k*D[2](T)(x, b) = h*(-T(x, b) + T3);

sol1 := simplify((pdsolve([pde, bc], T(x, y)) assuming (0 < a, 0 < b)));

#and the constant value are:  a = 250, b = 4, k = 2.091, T1 = -5, T2 = 0, h = 100, T3 = 1000

#the plot range is 0<x<250, 0<y<4.
params:={a = 250, b = 4, k = 2.091, T1 = -5, T2 = 0, h = 100, T3 = 1000};
sol1n:=eval(sol1,params union {infinity=100}); # Choosing infinity =100 (!)

plot3d(rhs(sol1n),x=0..250,y=0..4);

### An animation in infinity (an attempt at justifying infinity = 100):
plots:-animate(plot3d,[eval(rhs(sol1),params),x=0..250,y=0..4],infinity=10..100, frames=9);

The animation in the parameter 'infinity':

Notice the difference between subs and eval:
 

subs(x=0,exp(x));  #exp(0)
%; # 1
eval(exp(x),x=0); # 1

That is by design.

Using assumptions helps:
 

restart;

expr:=a^2*(2*a^2*p^3-a^2*((p^2+1)^2*(a^2-1))^(1/2)-((p^2+1)^2*(a^2-1))^(1/2)*p^2+2*p*a^2-2*p^3+((p^2+1)^2*(a^2-1))^(1/2)-2*p)/((p^2+1)^2*(a^2-1))^(1/2)/(p^3-((p^2+1)^2*(a^2-1))^(1/2)+p)/(a^2-p^2-1);

expr2:=simplify(expr) assuming a>0,p::real;
res:=int(expr2,p);
diff(res,p)-expr2;
simplify(%) assuming a>0,p::real;
#####
idts:=indets(res,name);
select(type,idts,'indexed');
op~(0,%);
indets(%,`local`);
c:=op(%);
resC:=subs(c[5]=C,res);

You see that the result res contains an escaped local 'a' in an indexed form as a[5] besides your global a.

Had you used a different name e.g. b this escaped local would still have been 'a'..

Anyway, even this shouldn't happen.

Clearly a bug.

restart;
ode:=(y(x)-x*diff(y(x),x))/(y(x)^2+diff(y(x),x))=(y(x)-x*diff(y(x),x))/(1+x^2*diff(y(x),x));
sol:=dsolve(ode);
### Isolating y'(x):
odes:=solve(ode,{diff(y(x),x)});
dsolve~([odes]); # Fine

 

Two exact solutions are easily found:

restart;
###
A1 := 8*Pi^3*R^2*n(x)^4*m+(2*Pi*sin((1/2)*x)*m*omega0*p+Pi*sin((1/2)*x)*m*omega0+3*Pi^2*(diff(n(x), x, x)))*n(x)^3+(-2*sin((1/2)*x)^2*m^2*omega0^2*p^2+2*cos((1/2)*x)^2*m^2*omega0^2*p^2-2*sin((1/2)*x)^2*m^2*omega0^2*p+2*cos((1/2)*x)^2*m^2*omega0^2*p)*n(x)^2+(-4*(diff(n(x), x$2))*sin((1/2)*x)^2*m^2*omega0^2*p^2-8*sin((1/2)*x)*(diff(n(x), x))*cos((1/2)*x)*m^2*omega0^2*p^2-4*(diff(n(x), x$2))*sin((1/2)*x)^2*m^2*omega0^2*p-8*sin((1/2)*x)*(diff(n(x), x))*cos((1/2)*x)*m^2*omega0^2*p)*n(x)+8*sin((1/2)*x)^2*(diff(n(x), x))^2*m^2*omega0^2*p^2+8*sin((1/2)*x)^2*(diff(n(x), x))^2*m^2*omega0^2*p;

###
R := 1; m := 1; p := 10; omega0 := 1000;
###
bc:=n(0) = n(4*Pi), D(n)(0) = D(n)(4*Pi);
###
### Ansatz: n(x) = A*sin(x/2)
eval(A1,n(x)=A*sin(x/2));
factor(%);
Asol:=solve(%/A^2,A);
### Thus we have the two exact solutions:
sol1,sol2:=Asol[1]*sin(x/2),Asol[2]*sin(x/2);
odetest(n(x)=sol1,A1); # 0
odetest(n(x)=sol2,A1); # 0
### bc is clearly satisfied for both solutions. 
### Plot:
plot([sol1,sol2],x=0..4*Pi);

You can rather easily get a larger image by using the following code, where the appropriate size depends on your monitor's screensize:

DGR:=DrawGraph(GH, stylesheet = "legacy"):
plots:-display(DGR,scaling=unconstrained,size=[1800,default]);

You can also replace 'default' by an appropriate number (e.g. 800).

You should certainly also look into the many options for DrawGraph.

The following copy has luckily been shrunk by MaplePrimes:

You can do like this:

restart;
Digits:=15:
sigma := 1 + x^2/2;
k := 0.1;
Q := 0.5516;
### Use unapply:
lambda1 := unapply(-3*(7*k^3*sigma^3 + 32*Q*k^2*sigma^2 - 11*k^2*sigma^3 + 54*Q^2*k*sigma - 44*Q*k*sigma^2 + 11*k*sigma^3 + 36*Q^3 - 54*Q^2*sigma + 32*Q*sigma^2 - 7*sigma^3)/(20*sigma^4),x);
###
N:=6/10/10^(-6);
M := Matrix(N + 1, 2, datatype = float[8]);
###
for i from 1 to N+1 do x:=0.6*(i-1)/N; M[i]:= <lambda1(x)| x > end do:
M;
M[-1]; # the last row

But this is much faster:
 

restart;
Digits:=15:
sigma := 1 + x^2/2;
k := 0.1;
Q := 0.5516;
lambda1 := unapply(-3*(7*k^3*sigma^3 + 32*Q*k^2*sigma^2 - 11*k^2*sigma^3 + 54*Q^2*k*sigma - 44*Q*k*sigma^2 + 11*k*sigma^3 + 36*Q^3 - 54*Q^2*sigma + 32*Q*sigma^2 - 7*sigma^3)/(20*sigma^4),x);
N:=6/10/10^(-6);
V:=Vector(N+1,i-> (i-1)*0.6/N,datatype=float[8]);
W:=evalhf(map(lambda1,V));
M:=<W|V>;

 

The following is simpler and works in Maple 12 too:
 

p:=piecewise(`&vartheta;`(tau)>=ϑl,1,0);
eq := diff(`&vartheta;`(tau), tau, tau)+6.666666666*sin(`&vartheta;`(tau))+66.66666666*cos(`&vartheta;`(tau))^2*p*(`&vartheta;`(tau)-.7227342478)+66.66666666*sin(`&vartheta;`(tau))^2*p*(`&vartheta;`(tau)-.7227342478);
`&vartheta;l` := .7227342478;
ic := `&vartheta;`(0) = 0, (D(`&vartheta;`))(0) = 8;
ld := dsolve([eq, ic], numeric, range = 0 .. 5);
odeplot(ld,[tau,p],refine=1,thickness=2, color=blue);

Assuming that the answer to my question about gln is that that is the vector given, then you can just eliminate the sine function.

Using solve and the abbreviations f and s for phi and sin we get:
 

restart;
sys:={f-3/2*s=0,-f^2+3/2*sqrt(1-s^2)=0};
sol:=solve(sys,[f,s]);
Sol:=allvalues(sol);
events:=[[Sol[1,1,1],halt],[Sol[2,1,1],halt]];
eval(sys,events[1,1]);
solve(%,s);
eval(sys,events[2,1]);
solve(%,s);

The events are then

events := [[f = 1/2*sqrt(2*sqrt(10) - 2), halt], [f = -1/2*sqrt(2*sqrt(10) - 2), halt]]

and as shown if either of these events occur then both of the equations in sys are satisfied and vice versa.

I agree that this should have been handled better by odetest.

Looking at the code for `tools/map` gave me the idea of using assumptions.
The following returns 0 in a very short time:

odetest(sol,ode) assuming real;

Another observation:

sol1:=dsolve(ode);
simplify(rhs(sol1)-rhs(sol)) assuming real;
eval(%,_C1 = C[1]);   # 0.

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