Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@tomleslie It should be pointed out that the numerical solution for U(x) found by using SYSxt4 is not differentiable at x=2 (still using x0=1). That is, if we take the numerical solution to be an approximation of U(x)=0 for x>2. Otherwise one must say that the solution is wrong for x>2.

I was somewhat surprised that there weren't any error messages from the numeric solver as there was when considering the original system and using rkf45_dae.

This system was discussed earlier in a related question by patient and he gave a symbolic solution.
http://www.mapleprimes.com/questions/204142-Error-in-DEtoolsconvertsys-Unable

restart;
SYSxt4:={diff(S(x), x) = -(-2*sqrt(S(x))*U(x)+2*S(x))/x, diff(U(x), x) = -U(x)*(-4*U(x)^2*sqrt(S(x))+8*U(x)*S(x)+3*S(x)^(3/2)*x^2-4*S(x)^(3/2))/(S(x)*x*(8*U(x)+x^2*sqrt(S(x))-4*sqrt(S(x))))};
#The symbolic solution is
SOL:={S(x) = exp(-(1/4)*x^2), U(x) = -(1/4*(x^2-4))*exp(-(1/8)*x^2)};
#Check:
odetest(SOL,SYSxt4);
simplify(%) assuming real;
ICS:={S(1) = exp(-1/4), U(1) = (3/4)*exp(-1/8), (D(S))(1) = -(1/2)*exp(-1/4)};
subs(x=1,ICS[1..2],SOL);
#The solution SOL of SYSxt4 and ICS[1..2] is certainly the only one in the open interval (0,2).
At x=2, U(x)=0, so the numerator and denominator of the expression for diff(U(x),x) in SYSxt4 are both zero, which clearly would be a problem for the numerical solver.
The exact solution has limiting derivative at x=2 equal to -exp(-1/2), i.e. not zero.

How about making use of the examples in the help page for dsolve:

?dsolve

@momoklala I guess we can just keep going on forever with changing values for the parameters.

Basically you must study and try the different options for dsolve/numeric/bvp.

Solving bvp problems can be challenging.

So how do you propose to get rid of diff(S(x),x,x)?

If somebody is going to do anything on this, he will surely like code that can be copied and pasted!

@momoklala With (now S=0.5  !!!) you just need to add the option initmesh=512 to the dsolve command that I gave earlier.

@momoklala The image you provided just proves to me that somehow this error was produced.
But I cannot even see the whole worksheet.

To get anywhere from here we shall need the actual worksheet, not an image of it.
To upload a worksheet (i.e.. the mw-file itself) use the thick green arrow in the editor to the right.

@Alejandro Jakubi Thanks for digging deeper into this. The difference between pp and qq is illuminating.

I submitted a bug report about convertsys giving the example diff(x(t),t,t) + diff(a(t),t,t) = 0.

Although surely the procedure was intended primarily for dsolve/numeric, it is not mentioned in the documentation that known functions can contain no derivatives.
I ran into this problem several years ago as mentioned above, but cannot remember if I filed a bug report at that time.

@StefKi Your system is linear in the second derivatives, so it ought to be simple to solve for those (when possible).
convertsys has to do that.

I tried doing it myself:

res:=solve(Sys,diff(GenKo,t,t)):
nops(res); #7 as expected
length(res); #Huge so do indeed use colons
FirstOrdSys := convertsys(res, [], GenKo, t, Y, YP ): #working on res instead of Sys
whattype(FirstOrdSys); #list
nops(FirstOrdSys); #4 as expected
length(FirstOrdSys); #Huge
length(FirstOrdSys[1]); #The bulk of the length
FirstOrdSys[2]; #The 'translation'
FirstOrdSys[3..4];

The total time used was less than 2 minutes.

Finally, your system doesn't contain derivatives of known functions (i.e. functions that are not dependent variables). Therefore the handling (or lack thereof) of such systems is irrelevant for your system.

@momoklala I have no problem with S:=-0.5.  I get

[-1.18614066271840, -1.17788387531395, -1.16926687825642, -1.16025994699062, -1.15082789257908, -1.13052356717881]

I use Maple 2015. I shall test it in Maple 18 later.
"Later": Well, it gives me the exact same results without any problem.

So what are you doing?

Several years ago I wrote a KonvSys procedure tailored for teaching purposes. I simply built it on top of `DEtools/convertsys`. Looking at it now I see that I ran into the very same problem raised here.
I solved it at that time by using the following modification, which still works:

STR:=convert(eval(`DEtools/convertsys`),string);
STR1:=convert('indets(ndeqns,'specfunc(anything,diff)')',string);
STR2:=convert('select(has,indets(ndeqns,'specfunc(anything,diff)'),u)',string);
konvsys:=parse(StringTools:-Substitute(STR,STR1,STR2));

Now just use konvsys instead of `DEtools/convertsys`.

@Alejandro Jakubi I believe you are wrong. Changing the local to something else like aa doesn't change the behavior. I maintain that the difference in behavior between using 'a' or 'p' is due to the ordering I mentioned.

Consider the following extremely stripped version of the procedure `DEtools/convertsys`. The problems appear early so it is not very long:

restart;
pc:=proc(deqns,nu)
local diffset,a,diffmax,idiff,j,ndeqns;
if type(deqns,`=`) then ndeqns:={deqns} elif type(deqns,'list') then ndeqns:={op(deqns)} else ndeqns:=deqns end if;
diffset:=indets(ndeqns,'specfunc(anything,diff)');
for a in nu do diffmax[a]:=0 end do;
printf("%a\n",Diffmax[a]=diffmax[a]);
printf("%a\n",convert(diffset,list));
for idiff in diffset do
  j:=0;
  printf("%a %a\n",'j'=j,idiff1 =idiff);
  while type(idiff,'specfunc(anything,diff)') and nops(idiff)=2 do
    j:=j+1;
    idiff:=op(1,idiff);
    printf("%a %a\n",'j' =j, idiffwhile = idiff)
  end do;
  printf("%a %a\n",idiff2=idiff,'op'(0,idiff)=op(0,idiff));
  diffmax[op(0,idiff)]:=max(j,diffmax[op(0,idiff)]);
  printf("%a\n",Diffmax[op(0,idiff)]=diffmax[op(0,idiff)])
end do;  
end proc:
##The two test equations:
test:=diff(x(t),t,t)+diff(a(t),t,t)=0;
testp:=diff(x(t),t,t)+diff(p(t),t,t)=0;
#Testing
pc(test,[x]);
pc(testp,[x]);
pc(test,[a]);
pc(testp,[p]);
##
## I'm convinced from this that the problem is as simple as exhibited right here in procedure pp:
pp:=proc(i,j) local dm; dm[a]:=max(i,dm[a]); dm[a]:=max(j,dm[a]) end proc;
pp(1,1);
pp(1,2);

#The problem surfaces in  pc(test,[x]) and not in the other 3 tests because of evaluation rules in procedures.

Notice that if you wrap max(j,diffmax[a]) in eval, then all 4 tests produce the error.



@StefKi In order to answer this we would have to see your system. You could upload a worksheet.

@StefKi What I meant by calling a(t) a known function as opposed to the unknown x(t) was that in your call to convertsys you are giving x(t) as the unknown function.

If a(t) is actually an unknown function and if you have several odes then just use convertsys on the system including a(t) as an unknown as in this simple example:

ode1:=diff(x(t),t,t)+diff(a(t),t,t)=0;
ode2:=diff(a(t),t,t)+x(t)*diff(a(t),t)=g(t);
DEtools[convertsys]({ode1,ode2},[],[x(t),a(t)],t,Y,YP);
#################
If your situation is not as above i.e. a(t) is not an unknown in that sense, then freezing and thawing might be what you need.
Again a simple example:

restart;
ode1:=diff(x(t),t,t)+diff(a(t),t,t)=0;
ode2:=diff(y(t),t,t)+x(t)=diff(b(t),t)+g(t);
vars:=[x(t),y(t)]; #The unknowns
AB:=remove(has,indets({ode1,ode2},specfunc(anything,diff)),vars);
ABeqs:=AB=~freeze~(AB);
thaw(ABeqs); #Just testing the thawing.
DEtools[convertsys](subs(ABeqs,{ode1,ode2}),[],vars,t,Y,YP);
thaw(%);
############
You might ask why some device like this freezing-thawing process is not part of the procedure convertsys.
My guess is that convertsys was primarily intended for numerical solution of odes, where conversion to a first order system is routinely done and where you will not encounter 'known' functions that are actually not known as is the case with your a(t).


@momoklala Try this:

restart;
Digits:=15:
Eq1 := diff(f(eta), eta, eta, eta)+f(eta)*diff(f(eta), eta, eta)-diff(f(eta), eta)^2 = 0;
blt := 10;
Eq2 := diff(theta(eta), eta, eta)/Pr+f(eta)*diff(theta(eta), eta) = 0;
bcs1 := f(0) = 0, D(f)(0) = 1, D(f)(blt) = 0;
bcs2 := D(theta)(0) = -N*(1+theta(0)), theta(blt) = 0;
q:=proc(n) local res;
  res:=dsolve(eval({Eq1, Eq2, bcs1, bcs2}, {Pr = 7,N=n}), numeric, output = listprocedure,initmesh=2048);
  subs(res,theta(eta))(0)
end proc;
plot(q,0.1..1.89,labels=[N,theta(0)]);
plot(q,0.1..2,adaptive=false,numpoints=20,labels=[N,theta(0)]);


Why do you say those values are wrong? What is your reason?

Incidentally, there is a missing assignment symbol (:=) in the assignment to bcs2.
Also it is not a very good idea to set Digits to 5. The default 10 should be fine, but if you would like to change then you should rather raise Digits to 15. However, that is not going to get you the values you (for some reason) want.

First 117 118 119 120 121 122 123 Last Page 119 of 231