Preben Alsholm

13743 Reputation

22 Badges

20 years, 341 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Aakanksha I tried the code in Maple 9.5 and Maple 10. Both gave results that were quite different from the results from Maple 11 and Maple 18. So something must have changed between Maple 10 and Maple 11.
The plot in Maple 18 (including the holes which are due to small imaginary parts caused by roundoff):

@Aakanksha Which Maple version are you using. In Maple 18 I tried your code (with or without RealDomain) but with Digits at default (10) and it worked fine. There were a few holes in the graph, but using Axel Vogt's suggestion I also tried to plot Re(c) and that graph turned up with no holes.

On this computer I happen to have access also to the much earlier Maple 11. The code worked there too.

restart;
k:=0.99: #with(RealDomain):
m:=1: #Digits:=2:
x:=(Pi*csc(Pi*(k-m)))/(0.693*GAMMA(k)*GAMMA(m));
m1:=MeijerG([[-m],[1-m]],[[0,-m,-m],[k-m]],(-m*k)/snr);
m2:=MeijerG([[-k],[1-k]],[[0,-k,-k],[m-k]],(-m*k)/snr);
c:=x*((((m*k)/snr)^m)*m1-(((m*k)/snr)^k)*m2);
with(plots):
plot(c,snr=0..10);


@belief111 You should take a look at the output of solve(eq,u) in the case you have in your worksheet allvalue.mw, where K:=1. In fact you do that after trying allvalues. You see that solve already gives you all the values although it is using indexed RootOfs. So by doing allvalues(solve(eq,u)) you are applying allvalues to a sequence. That means that allvalues misunderstands the input. Indeed, the 2nd root will be seen as a RootOf-option (see the help) and the 3rd through the 6th roots will be seen as possible other options to allvalues and as such fails miserably.
You may try (still with K:=1)
allvalues~([solve(eq, u)]); #Elementwise application of allvalues
but that doesn't (and shouldn't) give you anything new, but is syntactically correct.
Try after  K:='K'; to execute
solve(eq,u);
you will notice that the output is a single nonindexed RootOf. That represents all the 6 roots and those you can get by using allvalues.



@Joe Riel I only got Iterator.hdb and Iterator.mla in Iterator/lib.
The error message from with(Iterator);  is

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

In Maple 17, where the help works I get the same error message.

The version I installed (or thought I installed) is 2.3.4.


@Joe Riel I got the most recent version of the Iterator package from the Maple Application Center.
I have used an earlier version before. I'm using Windows Vista. In Maple 18 the libname seems quite right:
"C:\Program Files\Maple 18\lib", "C:\Users\alsholm\maple\toolbox\Iterator\lib", "."
That User directory is the place where Iterator.hbd and Iterator.mla are located.
However, ?Iterator doesn't bring up the help for Iterator, but something from combstruct. Also with(Iterator); just gives an error message.
In Maple 17, however, where libname still includes the same User directory, ?Iterator does bring up the help page for Iterator. But also here with(Iterator); doesn't work.
Any guess as to what I might be doing wrong?

@nm It is unsafe to use sum since in that case Maple will only execute the procedures a and b twice: once for n=0 and once for general n. The exceptional cases will not be caught by int. If you replace sum with add the procedures a and b will be called maxN + 1 times so that int will work on concrete values of n.
To see this you can insert a print statement in the procedure a:
a:=proc(n) print(n);
           int(f(x)*cos(n*freq*x),x=from_..to_) /denomC;
         end proc;


Your equation is missing two multiplication signs, they have been put in by DJJerome below.

@Chaimongkol I entered:
ode:=(4*Wi^2*t^3+4*t+4*sqrt(1-t^2)*De)*u(t)+(-8*sqrt(1-t^2)*De*t-4*L-4*De^2*t^3)*diff(u(t),t,t)=0;
infolevel[dsolve]:=5:
dsolve(ode);
And got a lot of printout. At one point it said:
-> Trying a Liouvillian solution using Kovacic's algorithm
   <- No Liouvillian solutions exists

Then it continued with something else, but didn't finish with anything before my patience ran out.




@baharm31 Saving w-values to a vector or matrix is easily done. Here is the simple example I used above.

restart;
pde:= diff(w(x, t),x,x)-diff(w(x,t),t)=0;
pi:=evalf(Pi): #Instead of just 3.14.
res:=pdsolve(pde,{w(x,0)=sin(x),w(0,t)=0,w(pi,t)=0},numeric);
res:-plot3d(t=0..1);
W:=subs(res:-value(output=listprocedure),w(x,t));
W(pi,0.1); #test
#To see what is going on m and n are fairly small.
#Increasing x- and t-values along rows and columns respectively.
t0:=0: t1:=1: a:=0: b:=pi: m:=10: n:=8:
h:=(b-a)/(n-1): k:=(t1-t0)/(m-1):
A:=Matrix(m,n,(i,j)->W(a+h*(j-1),t0+k*(i-1)));
A[..,3]; #The third column: fixed x
A[4,..]; #The fourth row: fixed t
plots:-matrixplot(A);


@baharm31 I tried to come up with an example with an integral where an exact solution to the problem is known.
Here is one which can be varied to some degree.
I determine ft such that u(x,t) = sin(x+t) is a solution to the pde and the initial and boundary conditions.
restart;
pde:=diff(u(x,t),t,t)=diff(u(x,t),x,x)+int(diff(u(x,t),x)*u(x,t),x=0..1/2)+ft;
eval(pde,u(x,t)=sin(x+t));
S:=solve(%,{ft});
pde:=subs(S,pde);
ibcs:={u(x,0)=sin(x),D[2](u)(x,0)=cos(x),u(0,t)=sin(t),u(1,t)=sin(1+t)};
res:=pdsolve(pde,ibcs,numeric,timestep=0.05,spacestep=0.05);
res:-animate(u(x,t)-sin(x+t),t=5);
res:-plot3d(t=5);
plot3d(sin(x+t),x=0..1,t=0..5);
res:-plot3d(u(x,t)-sin(x+t),t=5);
#The maximal discrepancy is about 4 percent as you can see.
Variations: As integrand you can try diff(u(x,t),x)^2, u(x,t)^2, u(x,t)^3, which work, but it is puzzling that u(x,t) gives an error:
Error, (in pdsolve/numeric/process_PDEs) selecting function must return true or false

The same error message comes with the integrand diff(u(x,t),x).
################################ Comment #####################
I tried a similar example using my own unsophisticated code which handles pdes of the form
diff(u(x,t),t) = f(x,t,u,ux,uxx)
where f is just a function of the 5 real variables (not functions) x,t,u,ux,uxx. The unamended code delivered a result on an f that included an integral of the function u(x,t)^2 w.r.t. x which was not at all ridiculous.
I conclude that the fact that pdsolve occasionally seems to handle integrals is unintended. Also there is no claim made in the help pages that it should do that.
The integral int(u(x,t)^2,x=0..c) is in my unamended code just replaced by c*u[j,n]^2. That clearly is not a very good substitute for the integral. It ought to be something like h*add(u[jj,n]^2,jj=0..j1), where j1 corresponds to the x-value c and where h is the spacestep.


@baharm31 I tried setting infolevel[int]:=5.
It produced a lot of printout, but ended claiming to have succeeded (to my great surprise).
Then I tried replacing the integral with 0 and there clearly v was just identically zero, which it ought to be.
After that I tried replacing the integral by its inegrand:
pde20:=evalindets(pde2,specfunc(anything,int),s->op(1,s));
res2:=pdsolve({pde1,pde20,pde3},bcs,numeric);
res2:-plot3d(w(x,t),t=0..tmax);
res2:-plot3d(v(x,t),t=0..tmax);
res2:-plot3d(u(x,t),t=0..tmax);
p:=res:-value(v);
p(xmax*.123456,tmax*.8765);
p2:=res2:-value(v);
p2(xmax*.123456,tmax*.8765);
#Clearly different, so maybe the integral is handled OK, but I'm still sceptical.
This time I haven't tried Dirac, but from what I said earlier that would have to be treated like 0.



@baharm31 After debugging pdsolve/numeric it appears that a check in the procedure `pdsolve/numeric/par_hyp` of the orders of the highest derivatives w.r.t. the space variable is performed.
It seems that it is necessary for the algorithm that the sum of the orders for the dependent variables (v and w) should be the same as the sum of the highest orders for each equation. In your case the first sum is 0+4 = 4 and the second sum is 2+4 = 6. Since 4<>6 you get the error message about not being close to standard form.
To make a small test of this try replacing diff(w(x,t),x,x) by w(x,t) in PDE1a. This will bring down the second sum to 0+4= 4, which makes it agree with the first:
res:=pdsolve({subs(diff(w(x,t),x,x)=w(x,t),PDE1a),PDE2a},bcs2,numeric);
res:-plot3d(w(x,t),t=0..tmax);
res:-plot3d(v(x,t),t=0..tmax);

#############Introducing a third dependent variable ###############
It appears that by introducing one more dependent variable the algorithm will work:
pde1:=diff(w(x,t),x,x)=u(x,t);
pde2:=subs(pde1,PDE1a);
pde3:=subs(pde1,PDE2a);
#Notice that the first sum is now (for u,v,w): 2+0+2 = 4.
#The second sum is (pde1,pde2,pde3): 2+0+2 = 4.
bcs:={v(x, 0) = 0, w(0, t) = 0, w(x, 0) = 0, D[1](w)(0, t) = 0, D[2](w)(x, 0) = 0, u(.57, t) = 0, D[1](u)(.57, t) = 0};
res:=pdsolve({pde1,pde2,pde3},bcs,numeric);
res:-plot3d(w(x,t),t=0..tmax);
res:-plot3d(v(x,t),t=0..tmax);
res:-plot3d(u(x,t),t=0..tmax);




@baharm31 bcs2 has the initial and boundary conditions for both of the pdes. Thus the syntax ought to be:

PDES := pdsolve({PDE1a,PDE2a}, bcs2,numeric);

However, Maple gives the error:

Error, (in pdsolve/numeric/par_hyp) input system is too far from a 'standard' form (see ?pdsolve,numeric for more detail)

Then when you go to that help page it says in the section with the heading 'Time-based Solver':

"The first mode of operation uses the default method, which is a centered implicit scheme, and is capable of finding solutions for higher order PDE or PDE systems. The PDE systems must be sufficiently close to a standard form for the method to find the numerical solution."

That doesn't really give us the promised "more detail".




@baharm31 You will have to get rid of the integral, but then you have an entirely different system.
To solve these two much simplified equations together numerically you have to replace v(t) by v(x,t) in both equations, which is easily done by subs:
PDE1a:=subs(v(t)=v(x,t),PDE1); #Similarly in PDE2
Secondly, you have to provide an initial condition for v(x,t) of the form v(x,0)=f(x), where f is known.
Then simply do
res:=pdsolve({PDE1a,PDE2a], bcs union {v(x,0)=f(x)},numeric); #Add options as you prefer

@baharm31 By without Dirac I thought that you meant to leave out the term having Diract in it. That term is the only one involving v.
So what do you mean by PDE2 without Dirac?

First 140 141 142 143 144 145 146 Last Page 142 of 231