Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

You are missing a value for 'u' for one thing.

Since you have an iterated integral you ought to see if you can get a value out of just e.g. h2(1.2345) with the value of u that you have. If that fails then forget about plotting.

@Marcel123 My general advice is: Do not ever assign to the variables used in odes, in your case x and t.
The variables are used in the solution as just variables and shouldn't be changed.

Use other letters, as in

X:=t->rhs(sol(t)[2]);   ## Capital X

@Z1493 From what I have seen so far my answer to the question "I only wanted to know if I have a PDE system with composite functions (as in the above case, f=f(x(t),t) ) how would I proceed to integrate it"
would be: Solve the pde system, then insert x=x(t).
But I think that what is needed now is an actual physical situation where this comes up.

@Z1493 Let me point out that
(1) The pdes are not the same.
(2) The solutions given by pdsolve are building blocks. See my previous reply with the infinite sum.

Question: Are you asking if given one system of (linear) pdes with coefficients depending on t and another where t is replaced by ln(x) in the cofficients only then the general solution for the two different systems will agree if we replace x by exp(t) in both?
Or expressed differently: Along the curve (x,t) = (exp(t),t) the solutions are the same?
That may be right, but is it interesting? Because of my point (2) above this is not disproved by your example (where you don't really test for that though: You need to replace x by exp(t)).

@Z1493 It is not clear to me what it is you want to solve.
So instead of your abstract formulation with A and B could you formulate the problem (i.e. the equations) concretely from the start.
You don't mean (I take it):

restart;
eq1:=D[2](f1)(exp(t), t) = f2(exp(t), t);
eq2:=D[1](f2)(exp(t), t) = -f1(exp(t), t);

and then hope to solve this for f1 and f2?

@Z1493 Your pde system could be rewritten as a second order pde in f1 (which I shall call f):
pde:=diff(f(x,t),x,t)+f(x,t)=0;
pdsolve(pde);
pdsolve(pde,build);
combine(%);
#By linearity any (convergent) infinite linear combination of solutions is again a solution:
S:=Sum(c(n)*exp(q(n)*x-t/q(n)),n=1..infinity);
pdetest(f(x,t)=S,pde); #Gives 0 as expected.
##Example
eval(S,{c(n)=1/n^4,q(n)=I*n}); #Picked c(n) for fast convergence
solR,solI:=op(evalc([Re,Im](%))); #If you have a complex solution then the real and imaginary parts are also solutions:
pdetest(f(x,t)=solR,pde); #Similar for solI
#For plotting purposes we can replace infinity by e.g. 10:
plots:-animate(plot,[eval([solR,solI],infinity=10),x=0..10],t=0..10,frames=100,size=[1000,400]);






We must expect that there are other nontrivial solutions corresponding to other values of L.
To find one more I found it easier to use a 2nd order ode satisfied by y(xi).

UU:=(isolate(sys[2],U(xi)));
eval(sys[1],UU);
ode:=collect((lhs-rhs)(%),diff,factor)=0;
odeL:=subs(y(xi)=L*y(xi),convert(ode,D));
resL:=dsolve({odeL,cond,D(y)(1)=1},numeric,approxsoln=[L=1,y(xi)=(1-xi)*(xi-0.8)]);
resL(.9); #Finds L = 1 approximately (as before)
UU; #(the right hand side can be plotted)
plots:-odeplot(resL,[[xi,y(xi)],[xi,rhs(UU)]]);
## Now looking for another value of L near 4:
resL2:=dsolve({odeL,cond,D(y)(1)=1},numeric,approxsoln=[L=4,y(xi)=(1-xi)*(xi-0.8)*(xi-0.9)]);
resL2(1); #L very close to 4
plots:-odeplot(resL2,[[xi,y(xi)],[xi,rhs(UU)]]);


@RomchikM In the code I gave in my answer above you will see the line
plots:-odeplot(res,[[xi,y(xi)],[xi,U(xi)]]);

This shows the graphs of y and U.

########
The problem you are facing is very similar to the following problem.

restart;
ode:=diff(y(x),x,x)+one*y(x)=0;
bcs:= y(0)=0, y(Pi)=0;

dsolve({eval(ode,one=1),bcs}); #infinitely many solutions (as is well known)
dsolve({eval(ode,one=1+1*10^(-8)),bcs}); #No nontrivial solution
res:=dsolve({eval(ode,one=1+1*10^(-8)),bcs},numeric); #Numerical attempt
plots:-odeplot(res,thickness=3);# same result 0
res2:=dsolve({eval(ode,one=1+1*10^(-8)),bcs},numeric,approxsoln=[y(x)=sin(x)]);
plots:-odeplot(res2,thickness=3); #Tiny, but it resembles a*sin(x) for a very small a
resL:=dsolve({ode,bcs,D(y)(0)=1},numeric);
resL(0); #'one' is very close to 1
plots:-odeplot(resL);



So I would suggest using the results I gave in my answer with the adjusted parameter L. Your input equations contain floats and are surely rounded, so why not?



 

@MTata94 Sure.
But you could also define the procedure q outside p, but then it shouldn't be declared local to p.
If you have a lot of subroutines used in p you may want to define them outside p and then wrap the whole thing in a module where the subroutines could be declared local. That is another story.

Thanks to all.
From the responses I concluded that submitting an SCR is justified. After all, SCR stands for "Software Change Request". It is not only (I hope) a euphemism for bug report.
The request is that dsolve/numeric with compile=true be made to work without a user with no particular technical computer skills having to do anything out of the ordinary.

I searched MaplePrimes for "compile" and found the following discussion

http://www.mapleprimes.com/questions/138276-Win7-And-Compiler

In that a different issue is raised about making Compiler:-Compile work in the Classical interface on 64 bit machines.

One of acer's workarounds redefines 'ssystem' to `dsolve/numeric/ssystem`.
That made me take a look at that procedure on my 64 bit machine with Windows 10 and Maple 2015.1:

I did as follows:

restart;
eval(`dsolve/numeric/ssystem`);
debug(`dsolve/numeric/ssystem`);
dsolve({diff(x(t),t)=x(t),x(0)=1},numeric,compile=true);

I got the following from the last command:

{--> enter dsolve/numeric/ssystem, args = cl -c -DX86_64_WINDOWS -D_M_X64 -we4028 -we4029 -I\"C:\\Program Files\\Maple 2015\\extern\\include\" /Oiyt -DMSVC -Gz -MD -Gy -nologo \"C:\\Users\\Bruger\\AppData\\Local\\Temp\\m2c98262.c\" -F\"oC:\\Users\\Bruger\\AppData\\Local\\Temp\\m2c98262.obj\", 1
                           "ode2.dll"
proc()  ...  end;
                            [1, ""]
<-- exit dsolve/numeric/ssystem (now in comp_proc) = [1, ]}
###
and then the error message reported earlier.

Actually, I don't understand any of it, but maybe somebody else does?
##
I located the file m2c98262.obj in the Temp directory given above. Opened it in Notebook. It is quite readable and starts with:
/**************************************************************
   External C module
   Code translation and wrapper code:
**************************************************************/




@Markiyan Hirnyk Yes, I should have added that on another computer, a 32 bit, on which I also installed Windows 10, I had no problems. This was using Maple 12 though. There Watcom is used.

@Carl Love Sorry, I'm using Maple 2015.1.

Your code for the loop is broken.
I would strongly advise you to use 1-D math input and worksheet interface whenever you want to do anything serious in Maple, i.e. if your sole purpose is not to make things look pretty.

Go to Tools/Options/Display/Input display and pick Maple Notation (that is 1D math notation).
Right after that continue in the Options dialogue box to the Interface tab. Go to Default format for new worksheets.
Choose Worksheet.
Finally press the Apply Globally button at the bottom.
These settings can always be undone if you don't like them.
Notice that changes only show up in new worksheets.

@Boby An iteration method that depends on successively integrating symbolically is in general going to break down very fast because of the inability to complete the integration.
Just take the case g(y) = y*ln(y).
The following will break down already with m=3 even if you replace YF(0) with the limit from the right.

restart;
g:=y->y*ln(y);
Y[0]:=exp(1); #Example chosen for relative simplicity
for m from 1 to 2 do
  YF:=unapply(Y[m-1],x);
  #Y[m]:=limit(YF(x),x=0,right)-int(t^2*(1/t-1/x)*g(YF(t)),t=0..x) assuming x>0,x<sqrt(6)
  Y[m]:=YF(0)-int(t^2*(1/t-1/x)*g(YF(t)),t=0..x) assuming x>0,x<sqrt(6)
end do;
evalf(Y[2]);
collect(%,ln);
#Compare with numerical solution:
ode:=diff(y(x),x,x)+2/x*diff(y(x),x)+y(x)*ln(y(x))=0;
res:=dsolve({ode,y(0)=exp(1),D(y)(0)=0},numeric);
plots:-odeplot(res,[x,y(x)],0..50);
plot(Y[2],x=0..sqrt(6));
plots:-odeplot(res,[x,y(x)-Y[2]],0..sqrt(6),thickness=3);


First 109 110 111 112 113 114 115 Last Page 111 of 231