## 12103 Reputation

16 years, 195 days

## piecewise...

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

## Forgot multiplication sign...

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

## You have to buy it!...

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...

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 get another error instead...

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).

Preben Alsholm

## Numeric solution requires numeric values...

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...

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
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

## Error message gone...

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

## Elementwise operations...

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

## Declare local; printlevel...

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

## Integral never zero...

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

## Use cat...

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

## Use the 'known' option...

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

## Plotting and animating...

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
﻿