Preben Alsholm

12103 Reputation

22 Badges

16 years, 195 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

Assuming that x0 < x1 < x2 < ...  and that your equations are not equations but just algebraic expressions containing the variable x, you can do:

L:=[seq([x<v[i+1],S[i]],i=0..29)];
p:=piecewise(op(map(op,L))):
plot(p,x=x[0]..x[30]);
 

Preben Alsholm

You forgot a multiplication sign after k. Maple interprets k(xxx) as an application of the function k to xxx.

So do

int((k*(1-x+.3*x*ln(x))+a)^2,x = s .. 1);

Preben Alsholm

GlobalOptimization does not come with Maple 13.

It is an extra you have to buy.

Since I don't have it either, I can reproduce the error message:

with(GlobalOptimization);
Error, invalid input: with expects its 1st argument, pname, to be of type {`module`, package}, but received GlobalOptimization

Preben Alsholm

 

This may do some of what you want:

restart;
h:=x->sqrt(x);
g:=x->piecewise(type(h(x), And(realcons,Not(infinity) ) ), h(x), 0);
g(-1);
g(2);
g(-infinity);
g(undefined);

##Notice that also a name will return 0:
g(a);
##If that is undesirable you can let the procedure return unevaluated in that case, as in

g:=x->piecewise(type(x,name),'procname(x)', type(h(x),And(realcons,Not(infinity))), h(x), 0);
g(v);
eval(%,v=9);
 

Preben Alsholm

I did this,

restart;
M[Star]:=2.*10^(30);
M[planet]:=2.*10^(27);
q:=(M[planet])/(M[Star]);
#Delta[p]:=max(H,abs(R-a));
Delta4[p]:=max(H^4,(R-a)^4);
Alpha:=unapply(piecewise(R<a, -q^2*6.67*10^(-11)*M[Star]/(2*R)*R^4/Delta4[p] ,R>a,q^2*6.67*10^(-11)*M[Star]/(2*R)*a^4/Delta4[p]) ,R,a);
PDE1:=diff(S(R,t),t)=3/R*diff(R^(1/2)*diff( S(R,t)*R^(1/2)-2*Alpha(R,a)*S(R,t)^(3/2)/(6.67*10^(-11)*M[Star])^(1/2),R),R);
IBC := {S(0, t) = 0, S(99, t) = 0, S(R, 0) = f(R)};
eval(PDE1,{H=1,a=2});
PDE2:=subs(Float(undefined)=0,eval(PDE1,{H=1,a=2}));
#This was done to avoid another error
#The Float(undefined) is caused by differentiating the piecewise defined function.
f:=R->(99-R)*R;
#Just an example of an f
pds1 := pdsolve(PDE2, IBC, numeric, S(R, t), time = t, range = 0 .. 99, spacestep = 0.1e-1, timestep = 0.25e-4);
pds1:-plot(t=0.1,numpoints=5);

 and got the error message

Error, (in pdsolve/numeric/plot) unable to compute solution for t>0.:
unable to store -70000.0000000000+23165.0897338412*I when datatype=float[8]

indicating that S(t) is getting negative (S(t)^(3/2) is in the equation).

Let me add that I have not analyzed your system.

Preben Alsholm
 

You wrote that you want to solve PDE1 for a and H unknown.

Since the pde must be solved numerically, a and H must be given numeric values.

Preben Alsholm

I made these 3 procedures for my own teaching purposes (PolarForm, `print/Exp`, `value/Exp`):

PolarForm:=proc(tal::{algebraic,list,set},e::name:='useexp')
local EXP;
if type(tal,{list,set}) then return map(procname,tal,e) end if;
if type(tal,specfunc(anything,polar))
   then
      if e='useExp' then op(1,tal)*Exp(I*op(2,tal)) else op(1,tal)*'exp'(I*op(2,tal)) end if
   else
procname(radnormal(polar(tal)),e)
end if
end proc:
`print/Exp`:=proc(z) 'e'^(z) end proc:
`value/Exp`:=proc() exp(args) end proc:
PolarForm(1+I*sqrt(3));
%;
PolarForm(1+I*sqrt(3),useExp);
value(%);
 

Preben Alsholm

Maybe the following changed version will help. At least it doesn't produce an error message.

I have used another Delta, Delta4, to avoid abs. Also I have made Alpha a function with unapply

restart;
M[Star]:=2.*10^(30);
M[planet]:=2.*10^(27);
q:=(M[planet])/(M[Star]);
#Delta[p]:=max(H,abs(R-a));
Delta4[p]:=max(H^4,(R-a)^4);
Alpha:=unapply(piecewise(R<a, -q^2*6.67*10^(-11)*M[Star]/(2*R)*R^4/Delta4[p] ,R>a,q^2*6.67*10^(-11)*M[Star]/(2*R)*a^4/Delta4[p]) ,R,a);
PDE1:=diff(S(R,t),t)=3/R*diff(R^(1/2)*diff( S(R,t)*R^(1/2)-2*Alpha(R,a)*S(R,t)^(3/2)/(6.67*10^(-11)*M[Star])^(1/2),R),R);

Preben Alsholm
 

Say your vector is

v:= <a(t),b(t),c(t)>;

then in Maple 13 you can differentiate elementwise like this

diff~(v. t):

where you should notice the tilde.

In earlier versions:

map(diff, v, t);

map is still useful, though.

Preben Alsholm

Your procedure with two changes, viz. the local declaration and contourplot replaced by its long name:

PlotMfield := proc (a::float, b::float, c::float, d::float, p::float, chi::float)
local g,G,f,Bz,Bphi,Br;
g := z -> c*(exp(-(1/2)*z^2/a^2)/a-exp(-(1/2)*z^2/b^2)/b)/sqrt(2*Pi) ;
G :=z -> int(g(zeta), zeta = -infinity .. z);
f  := r -> d*exp(-p*r);
Bz := chi*p*(p-1/r)*G(z)*f(r);
Bphi := p*g(z)*f(r);
Br := chi*p*g(z)*f(r);
plots:-contourplot(sqrt(Br^2+Bphi^2+Bz^2), r = 0 .. 15, z = (-1)*2.5 .. 2.5, grid = [15, 75], scaling = constrained);
end proc;
 

Your second problem can be solved by executing
printlevel:=2; 
(because of the nested loop)
or by giving an explicit print command as in

for i in La do
   for j in Lc do
     print(PlotMfield (i,1.1*i,j,1.0,0.4,1.0))
   end do;
end do;

Preben Alsholm
 

The integral is never zero, and you need to change 'pi' to 'Pi'.

You can try the following instead.

F:=(t,E)->Int(sqrt((1-(1-E^2)*sin(x)^2)/(1-(1-E^2)*t^2*sin(x)^2)),x=0..Pi/2);
plots:-contourplot(F(t,E),t=0..1,E=0..1);
 kk:=evalf(F(.3,.4));
plots:-implicitplot(F(t,E)=kk,t=0..1,E=0..1);
 

Preben Alsholm

You could use concatenation, i.e. cat in Maple, as in

cat("Figure ", j+1);

I couldn't test it on your procedure since several things are missing (and you seem to be solving for a constant in b:=[solve(t^3-a=0,a)]);

Preben Alsholm
 

You can use the 'known' option, as in the example below:

dsolve({f1(0)=0,diff(f1(x),x)=f1(x)+1},numeric,output=listprocedure);
F1:=subs(%,f1(x));
g:=x->2*F1(x/2);
g(1);
dsolve({f2(0)=1,diff(f2(x),x)=sin(f2(x))+g(x)},numeric,output=listprocedure,known=F1);
F2:=subs(%,f2(x));
plot(F2,0..3);
 

Preben Alsholm

You didn't assign values to the constants, but just wrote equations.

That is one reason why no plot turned up. The other is that you are trying to plot u(x,t) as a function of 2 variables using 'plot'.

I wouldn't assign values to the constants myself, but would use eval as seen below.

plot(eval(u(x, 2), {F = 1, c = 1, l = 1, m = 1, omega = 1}), x = 0 .. 20);
plots:-animate(plot,[eval(u(x, t), {F = 1, c = 1, l = 1, m = 1, omega = 1}), x = 0 .. 20],t=0..5);

One small detail: The line
WaveEqn1 := subs(WaveEqn = u(x, t)):
doesn't do anything. Your assignment of u has already done what you intend to do.

Preben Alsholm

It seems to me that your code works. The warning is just that. A warning, not an error. It is to be expected since the function xx(t-1) at the time is undefined for t > 2.

Using output=piecewise makes it unnecessary to rename, since there is no recursive calling. In fact you made the right hand side of the differential equation into a concrete function of t given by a formula.

So you could change your code to

restart;
xx := t -> 1  :
P := plot(xx, -1 .. 0):
xnum := dsolve({x(0) = xx(0), (D(x))(t) = -xx(t-1)},type = numeric, output = piecewise, range = 0 .. 1):
xx:=unapply(subs(op(xnum),x(t)),t):
P := P, plot(xx, 0 .. 1):
#Just to see that it works
plots:-display(P);
xnum := dsolve({x(1) = xx(1), (D(x))(t) = -xx(t-1)},type = numeric, output = piecewise, range = 1 .. 2):
xx := unapply(subs(op(xnum),x(t)),t):
#Just to see that it still works
P := P, plot(xx, 1 .. 2):
plots:-display(P);
 

Preben Alsholm

First 147 148 149 150 151 Page 149 of 151