Preben Alsholm

13471 Reputation

22 Badges

20 years, 298 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Rouben Rostamian  And the logarithmic plot:

plots:-logplot(abs(myx1(t) - myx2(t)), t=0..4*Pi);


@GeorgeReynolds 
1. Consider the following simple example:
restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x,x);
pdsolve(pde,{u(x,0)=x,u(0,t)=0,D[1](u)(1,t)=1},numeric); ## OK
pdeD:=convert(pde,D);
pdsolve(pdeD,{u(x,0)=x,u(0,t)=0,D[1](u)(1,t)=1},numeric);

Error, (in pdsolve/numeric/process_PDEs) number of dependent variables and number of PDE must be the same

This is clearly a limitation in the implementation of pdsolve/numeric. It is stated explicitly in the help:
?pdsolve,numeric
"The PDE system PDEsys must be specified in diff notation."

There is no complaint from the symbolic solver:
pdsolve(pdeD); #Separates variables
pdsolve({pdeD,u(x,0)=x,u(0,t)=0,D[1](u)(1,t)=1}); #No error, but no result (as I expected)

2. No, surely not. I'm pretty certain that it uses discretization in space as well as in time.



@GeorgeReynolds You cannot use numerical solution unless the parameters are given numerical values. You have 4 parameters with no values given. One is I which you surely don't mean to be the imaginary unit that Maple takes it to be.
There are several solutions to this problem. The simplest (to me) is to rename your I to II or whatever you like.
There are 3 more:

indets(PDE1,name) minus {xi,vartheta};

Output:      {E, a, omega}

########### After that you will run into other problems. To avoid one of them correct your line
convert(PDE1, diff);
## to
PDE1 := convert(PDE1, diff);
## Next problem: Maple's numerical solver is time based. Let us assume that vartheta is your time. Thus it starts from initial values given on the xi interval xi=0..1 (in your case) and proceeds in time from there satisfying whatever boundary conditions are given at xi=0 and xi=1. You don't have any initial values, only boundary conditions.

##Example, where I have replaced I with II:
PDE1 := convert(PDE1, diff);
params:=indets(PDE1,name) minus {xi,vartheta};
parvals:=params=~1; #Setting all 4 parameters to 1
eval(PDE1, parvals); #Just to have a look
#ic:={V(xi,0)=xi,D[2](V)(xi,0)=0}; #An example
ic:={V(xi,0)=xi^2*(1-xi),D[2](V)(xi,0)=0}; #Maybe a better example
pds := pdsolve(eval(PDE1, parvals), bc union ic,numeric);
pds:-plot(vartheta=1);

 

@GeorgeReynolds To use Worksheet interface and 1D math input (in the session or globally) go to (in Windows):


Tools/Options/Display/Input display. Choose Maple Notation.  (i.e. 1D math).
Tools/Options/Interface/Default format for new workswheets. Choose Worksheet.

Finally press either Apply to Session or Apply Globally.

The changes only affects new worksheets.
There is no danger in choosing the latter since the change can always be reverted.  

What I see is (for the initial value problem, as an example):

y(x) = JacobiSN((1/2)*sqrt(2)*x+InverseJacobiSN(1, I), I)

and after a % I get

y(x) = JacobiSN((1/2)*sqrt(2)*x+EllipticK(I), I)

How did you get those sn's in the first place? Is there in fact something missing in your worksheet?

@Markiyan Hirnyk Since you need a restart in between, it is difficult to compare two exact result that are so big:

But just try

restart;
int(......);
evalf(%);

restart;
int(......);
evalf(%);

There you surely don't have numerical integration, so the cause is roundoff.

@Konstantin@ Yes, if int cannot do the integration (as is the case here), then it uses numerical integration.

Try

int(exp(sin(x)+x+exp(x)),x=0..1);

@Bendesarts The code I used:

restart;
K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);
omega[sw]:=beta/(1-beta)*omega[s];
for i to 4
do
r[i]:=sqrt((u[i](t)/(L/2))^2+(v[i](t)/H)^2):
omega[i]:=omega[st]/(1+exp(b*v[i](t)))+omega[sw]/(1+exp(-b*v[i](t))):
Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)+omega[i]*(L/2)/H*v[i](t):
Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)+omega[i]*(L/2)/H*v[i](t)+(K.<seq(v[i](t),i=1..4)>)[i]:
EqSys[i]:=[Equ[i],Eqv[i]]:
end do;
paramsGeo:=L=0.015,H=0.015,beta=0.5,Vf=0.3;
omegaS:=eval(Pi*Vf/L, [paramsGeo]);
paramsCycle:=omega[s]=omegaS,Au=1,Av=1,b=100;
params:=paramsGeo,paramsCycle;

sys:=map(op,eval({seq(EqSys[i],i=1..4)},{params}));
#ic:={u[1](0)=0.8, v[1](0)=0,u[2](0)=0.8, v[2](0)=0,u[3](0)=0.8, v[3](0)=0,u[4](0)=0.8, v[4](0)=0};
ic:={u[1](0)=0.8, v[1](0)=0.1,u[2](0)=0.8, v[2](0)=0.1,u[3](0)=0.8, v[3](0)=0.1,u[4](0)=0.8, v[4](0)=0.1};
res:=dsolve(eval(sys,omega[st]=50) union ic,numeric);
plots:-odeplot(res,[u[1](t),v[1](t)],0..1);
plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..1,thickness=2);
#The last plot:


I don't see any text in the body of the question. How come Carl can see it?

@GambiaMan What is it you want plotted: Only the points corresponding to the times in your sample? If that is the case your option numpoints=2000 confuses me.

@Bendesarts Change v to

vt:=<seq(a[i](t),i=1..4)>;

Without changing v:

v:=Vector(4,symbol=a);
vt:=apply~(v,t);

Comment: Why not just v(t)? Why this funny apply, which is otherwise almost never used?
Because () as well as [] are used these days for getting elements from vectors and matrices.
Thus v(1) returns a[1], and v(t) is understood as the t'th element of v, which doesn't make sense if t is not an integer.
There are differences between () and [], though. Look for programmer indexing contra mathematical indexing.
See
?Indexing Arrays, Matrices, and Vectors

## Making a function V which returns vectors:
V:=t->apply~(v,t);
V(s);
##Or the slightly strange looking version of the latter:
V:=unapply(apply~(v,t),t);
V(6);

@red1eco With the values p=50, q=0.01, W=0.1, I'm not surprised to find no visible change in the graph of x when tau is varied from 0 to 0.2.
In fact I made an animation in the parameter tau from 0 to 0.2. This animation wasn't very animated. For the fun of it I tried tau=0..20. That was much more exciting.
##At the bottom I have looked at differences instead. For that purpose I added the option output=listprocedure in dsolve below.
Note: In the first version I used the name p for the procedure I now call P (capital P). Since p is one of the parameters that had to be changed.


In your last worksheet to show the different graphs in the same picture you just need to save each plot in a variable right after it is produced, like this (e.g.) for the first one

plots:-odeplot(res, [t, x(t)], 0 .. 9, linestyle = dot);  p1 := %:

If you then have made the plots p1,p2,p3,p4,p5 you can display them together using display from the plots package:

plots:-display(p1,p2,p3,p4,p5);

You will be disappointed, however, as I warned you above: they fall almost on top of each other.

##########
The animation and the parameters option:
restart;
sys2 := [diff(x(t), t) = r*x(t)*(K-x(t)-a*y(t))/K, diff(y(t), t) = s*y(t)*(L-y(t)+b*x(t))/L-q*y(t)*E(t), diff(E(t), t) = W*((p-tau)*q*y(t)-c)*E(t)];
params := convert(indets(sys2, name) minus {t}, list); #LIST
paramvalsFix := {K = 2000, L = 1200, W = .1, a = 1.2, b = .2, c = 100, p = 50, q = 0.1e-1, r = 2, s = 1.2};
#tau left out above since we want to vary that
res:=dsolve({sys2[], E(0) = 75, x(0) = 451, y(0) = 1290}, numeric,parameters=params,output=listprocedure);
res(parameters=[op(paramvalsFix),tau=0]); #Setting parameter values
plots:-odeplot(res, [t, x(t)], 0 .. 9, linestyle = dot); #Example plot
## Now setting up a procedure to do the parameter setting and to produce the plot.
##The procedure allows optional arguments and those will be passed to odeplot (the _rest inside P).
P:=proc(tt) res(parameters=[op(paramvalsFix),tau=tt]); plots:-odeplot(res, [t, x(t)], 0 .. 9,_rest) end proc;
P(0); #Test with tau=0
P(0,color=blue); #Test with tau=0 and using the color blue
##Animation using p:
plots:-animate(P,[tau,caption="Superply dull"],tau=0..0.2);
plots:-animate(P,[tau,caption="More exciting and blue",color=blue],tau=0..20);
##The following indicates how you can display several graphs of x using the procedure P:
plots:-display(P(0,linestyle=dot),P(5,linestyle=dash),P(20,color=blue));
########DIFFERENCES:
X:=subs(res,x(t)); #Saving the procedure for computing x(t) in X (capital X).
N:=200: #Number of t-points
T:=<seq(9.0*i/N,i=0..N)>; #t-point vector
T[-1]; #The last point
##Now a procedure that produces a corresponding vector of x-values:
Q:=proc(tt) res(parameters=[op(paramvalsFix),tau=tt]); X~(T) end proc;
plot(T,Q(0)-Q(0.05),caption="The difference between x values for tau=0 and tau=0.05");





Maybe replacing sqrt(Drag*g/m) by sqrt(g/m)*drag2. The first appearance of Drag should then be replaced by drag2^2.

I have no reasonable way of testing it, since I don't have your data.

Imaginary values could appear this way:

ln(cosh(sqrt(-65.)));
                        -1.576130914+3.141592654*I

@red1eco In order to continue you can then do

paramvals := {K = 2000, L = 1200, W = .1, a = 1.2, b = .2, c = 100, p = 50, q = 0.1e-1, r = 2, s = 1.2, tau = .1};
res:=dsolve(eval({sys2[],x(0)=1.2,y(0)=0.5,E(0)=.1},paramvals),numeric);
plots:-odeplot(res,[t,x(t)],0..1);
plots:-odeplot(res,[[t,x(t)],[t,y(t)],[t,50*E(t)]],0..1,thickness=3,legend=[x,y,E]);


The same initial values have been used.

(I use my real name as user name. I retired from the Technical University of Denmark in 2011.)

@acer To compare the 3 methods I tried with an ode having a symbolic solution:

ode:=diff(x(t), t, t)+diff(x(t), t)+x(t) = 0;

Using the same initial values x(0)=-1,D(x)(0)=-0.8 and evaluating the second derivative at t = 4*Pi/sqrt(3).

The exact solution is easily found by

res:=dsolve({ode,x(0)=-1,D(x)(0)=-0.8});

#to be res:=x(t) = -13/15*sqrt(3)*exp(-1/2*t)*sin(1/2*sqrt(3)*t)-exp(-1/2*t)*cos(1/2*sqrt(3)*t)
##The result for x '' (t) at  t = 4*Pi/sqrt(3) is
rE:=eval(diff(rhs(res),t,t),t=4*Pi/sqrt(3));
r4:=evalf(rE):

The differences r1-r4, r2-r4, r3-r4 for the 3 methods (with Digits=10) were:

-0.301619e-5, -8.60919101217106*10^(-8), 9.21698134481730*10^(-12)

I also tried Digits = 15. That was disastrous for the plot in method 1 and the results were:
-.289843880257555, -8.61594650208852*10^(-8), -5.83379178298316*10^(-11)

It has been pointed out (coming indirectly from Allan Wittkopf, I think) that in general it is best to let dsolve/numeric do as much of the computation as possible.
This certainly confirms that.




First 91 92 93 94 95 96 97 Last Page 93 of 225