Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

The initial conditions are allowed to be discontinuous at a point:

eq:=diff(pA(z,t),t)=-diff(pA(z,t),z);

ibc:= [pA(z,0) = 0, pA(0,t) = 1];

res:=pdsolve(eq,ibc,numeric,time=t,range=0..1,spacestep=1e-3);

res:-animate(z=0..1,t=0..1,numpoints=300);
res:-plot3d(z=0..1,t=0..1,axes=boxed);

This simple equation is easily solved:

pdsolve(eq);
with result pA(z, t) = _F1(-z+t). The conditions ibc then imply

pA(z,t) = piecewise(t<z,0,1) (i.e. Heaviside(t-z))

plot3d(piecewise(t<z,0,1),z=0..1,t=0..1,axes=boxed);




The initial conditions are allowed to be discontinuous at a point:

eq:=diff(pA(z,t),t)=-diff(pA(z,t),z);

ibc:= [pA(z,0) = 0, pA(0,t) = 1];

res:=pdsolve(eq,ibc,numeric,time=t,range=0..1,spacestep=1e-3);

res:-animate(z=0..1,t=0..1,numpoints=300);
res:-plot3d(z=0..1,t=0..1,axes=boxed);

This simple equation is easily solved:

pdsolve(eq);
with result pA(z, t) = _F1(-z+t). The conditions ibc then imply

pA(z,t) = piecewise(t<z,0,1) (i.e. Heaviside(t-z))

plot3d(piecewise(t<z,0,1),z=0..1,t=0..1,axes=boxed);




'dsolve' should be 'pdsolve' as in the original question. I must have copied it wrong. I have corrected it above.

'dsolve' should be 'pdsolve' as in the original question. I must have copied it wrong. I have corrected it above.

The definition of B does not agree with ode[1], and in the definition of srPP you probably want srP instead of hrP.

Do you have any reason whatsoever why s1 and s2 (as defined last!) has anything to do with the difference between the unknown exact solution and the one found numerically by the method rk45 with abserr=10^(-10),relerr=10^(-10)??

@kinesimario How is Y (capital Y) defined?

@kinesimario How is Y (capital Y) defined?

@kinesimario As you write it there is nothing wrong once you have corrected ic2 as explained by Jarekk.

But why use a stiff solver?

deq2 := diff(x(t),t,t) = -x(t)/(x(t)^2+y(t)^2)^(3/2), diff(y(t), t,t) = -y(t)/(x(t)^2+y(t)^2)^(3/2);

ic2 := x(0) = 1, D(x)(0) = 0, y(0) = 0, D(y)(0) = 1;

rosen := dsolve({ic2, deq2}, type = numeric, range = 0 .. 2*Pi);

plots:-odeplot(rosen,[x(t),y(t)],0..2*Pi);

You may find some fun in trying the following:

rosen1 := dsolve({ic2, deq2}, numeric);
plots:-odeplot(rosen1,[x(t),y(t)],0..2*Pi,scaling=constrained);
?dsolve,interactive,numeric
rosen1(initial);
rosen1(initial=[1.5,0,0,1]);
plots:-odeplot(rosen1,[x(t),y(t)],0..33,scaling=constrained);
rosen1(initial=[2,0,0,1]);
plots:-odeplot(rosen1,[x(t),y(t)],-500..500,scaling=constrained);

@kinesimario As you write it there is nothing wrong once you have corrected ic2 as explained by Jarekk.

But why use a stiff solver?

deq2 := diff(x(t),t,t) = -x(t)/(x(t)^2+y(t)^2)^(3/2), diff(y(t), t,t) = -y(t)/(x(t)^2+y(t)^2)^(3/2);

ic2 := x(0) = 1, D(x)(0) = 0, y(0) = 0, D(y)(0) = 1;

rosen := dsolve({ic2, deq2}, type = numeric, range = 0 .. 2*Pi);

plots:-odeplot(rosen,[x(t),y(t)],0..2*Pi);

You may find some fun in trying the following:

rosen1 := dsolve({ic2, deq2}, numeric);
plots:-odeplot(rosen1,[x(t),y(t)],0..2*Pi,scaling=constrained);
?dsolve,interactive,numeric
rosen1(initial);
rosen1(initial=[1.5,0,0,1]);
plots:-odeplot(rosen1,[x(t),y(t)],0..33,scaling=constrained);
rosen1(initial=[2,0,0,1]);
plots:-odeplot(rosen1,[x(t),y(t)],-500..500,scaling=constrained);

If u is actually a function of two variables x and t, this construction won't work, but you can do any of the following.

restart;
DE:=diff(x(t),t)=cos(x(t));
IC:=x(0)=0;
u:=(x,t)->x^2+sin(t);
s := dsolve( {DE,IC }, numeric, output = listprocedure );
X:=subs(s,x(t));
U:=t->u(X(t),t);
plot( U, 0..10 );
#A somewhat clumsy alternative:
U:=tt->subs(s,u(x(t)(tt),tt));
plot( U, 0..10 );

#If you use output=operator the clumsy version becomes less so:

restart;
DE:=diff(x(t),t)=cos(x(t));
IC:=x(0)=0;
u:=(x,t)->x^2+sin(t);
s := dsolve( {DE,IC }, numeric, output = operator );
X:=subs(s,x);
U:=t->u(X(t),t);
plot( U, 0..10 );
#Alternatively, and not so clumsy this time:
U:=t->subs(s,u(x(t),t));
plot( U, 0..10 );

If u is actually a function of two variables x and t, this construction won't work, but you can do any of the following.

restart;
DE:=diff(x(t),t)=cos(x(t));
IC:=x(0)=0;
u:=(x,t)->x^2+sin(t);
s := dsolve( {DE,IC }, numeric, output = listprocedure );
X:=subs(s,x(t));
U:=t->u(X(t),t);
plot( U, 0..10 );
#A somewhat clumsy alternative:
U:=tt->subs(s,u(x(t)(tt),tt));
plot( U, 0..10 );

#If you use output=operator the clumsy version becomes less so:

restart;
DE:=diff(x(t),t)=cos(x(t));
IC:=x(0)=0;
u:=(x,t)->x^2+sin(t);
s := dsolve( {DE,IC }, numeric, output = operator );
X:=subs(s,x);
U:=t->u(X(t),t);
plot( U, 0..10 );
#Alternatively, and not so clumsy this time:
U:=t->subs(s,u(x(t),t));
plot( U, 0..10 );

You are more likely to get good answers if you provide us with the details, e.g. in the form of an uploaded worksheet or the complete lines of code pasted into your posting.

I don't have that problem in Maple 15. In Maple 12, however, I do. But if you change the range to 0..2*Pi as in

odeplot(res,[e1r(t),e1i(t)],0..2*Pi,refine=1);

the error doesn't appear. Alternatively, change the range in the assignment to res like this:

res:=dsolve(s1 union s2,numeric, method =rkf45_dae,abserr=10^(-10),relerr=10^(-10),range=0..Pi,maxfun = 0,output=listprocedure);

I don't have that problem in Maple 15. In Maple 12, however, I do. But if you change the range to 0..2*Pi as in

odeplot(res,[e1r(t),e1i(t)],0..2*Pi,refine=1);

the error doesn't appear. Alternatively, change the range in the assignment to res like this:

res:=dsolve(s1 union s2,numeric, method =rkf45_dae,abserr=10^(-10),relerr=10^(-10),range=0..Pi,maxfun = 0,output=listprocedure);

First 202 203 204 205 206 207 208 Last Page 204 of 231