Preben Alsholm

MaplePrimes Activity


These are answers submitted by Preben Alsholm

To get Vector/Matrix/Array output from dsolve/numeric use output=Array(...).
Example:
 

restart;
sys:=diff(x(t),t)=y(t),diff(y(t),t)+x(t)=sin(t);
ics:=x(0)=0,y(0)=1;
sol:=dsolve({sys,ics});
b:=10: N:=64: h:=b/(N-1);
T:=Array(1..N,i->i*h,datatype=float);
res:=dsolve({sys,ics},numeric,output=T);
A:=res[2,1];
plot(A[..,1..2]); # Test
X:=A[..,2]; # X values as a vector
Y:=A[..,3]; # Y values as a vector
plot(T,Y); # Another test

 

For the boundary conditions you have Maple's dsolve/numeric/bvp has severe problems and that may be because the problem has no solution. For the conventional boundary values f(0)=0, D(f)(0)=0, D(f)(inf) = k (k some constant, e.g. 2) there is no problem:
 

restart;
ode:=diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta$2)=0
bcs1:=D(f)(0)=1,(D@@2)(f)(0)=0,(D@@2)(f)(6)=2; #Your boundary conditions
dsolve({ode,bcs1},numeric,initmesh=8192); # No success even with initmesh high
## Conventional boundary conditions:
bcs2:=f(0)=0,D(f)(0)=0,D(f)(6)=2;
res2:=dsolve({ode,bcs2},numeric);
plots:-odeplot(res2,[seq([eta,diff(f(eta),[eta$i])],i=0..2)]);

Since dsolve/numeric/bvp has problems I'm not at all surprised that your scheme does too. On the contrary that should be expected.
To examine your scheme you may try using the boundary conditions bcs2 instead. If your scheme doesn't work quite as then look for logical mistakes.

The difference between the 'good' and the 'bad' solves the corresponding homogeneous pde.
Both solutions are correct.
 

pdetest(good,eqn_sys); # OK
pdetest(bad,eqn_sys); # OK
##
solhom:=y(t,x)=op(combine(expand(rhs~(good-~bad))));
HOM:=D[1, 1](y)(t, x)+y(t, x) = 0;
pdetest(solhom,HOM); # OK

 

Although Kitonum has already answered your question perfectly fine I shall add mine too.

Your question may be related to your question of June 19 2018:
https://mapleprimes.com/questions/224941-Error-in-Dsolvenumericbvpconvertsys#answer250475 

Now you give us 3 odes with only 2 unknown functions. Two of the odes are linear with constant coefficients thus basically trivial:
 

ode1:=1.857142857*10^11*diff(u(x), x, x)=0;
ode2:=1.547619048*10^10*diff(w(x), x, x, x)-1.300000000*10^11*diff(w(x), x)=0;

In addition we have the homogeneous boundary conditions:
 

bcs:= u(0)=0, u(0.001)=0, w(0)=0, w(0.001)=0, D(w)(0)=0;

Using dsolve on that system:

dsolve({ode1,ode2,bcs});

gives you (not surprisingly) u(x)=0, w(x)=0.
But you have a third ode (call it ode3) that only involves w. That one doesn't have the trivial solution w(x) = 0 for all x:

eval(ode3,w=0);

gives you .3693149535 = 0.


 

Copying the examples from the CUDA help page I get at the given stage in Maple 2018.1 as well as in Maple 2017.3:

tCUDA := time[real](M1.M2);
Error, (in unknown) CUBLAS internal error

In Maple 2016.2, 2015.2, 18,  17, 16, and 15 there is no problem with the example section for CUDA on my machine. Those releases of Maple are all installed on the same computer.

I shall submit an SCR.

Although as I pointed out in a comment, x>a and x<>c evaluates to x>a the equivalent formulation
 

x>a and (x>c or x<c);

does not. It returns unevaluated.
And so does:
 

x>-infinity and x<infinity and (x>1/2 or x<1/2);

solve handles the latter just fine:

solve(x>-infinity and x<infinity and (x>1/2 or x<1/2),x);

and returns RealRange(Open(1/2), Open(infinity)), RealRange(Open(-infinity), Open(1/2)).

In general x<>c is not equivalent to (x<c or x>c) since the latter assumes an order where the first simply means 'x is different from c'.

As mehdi jafari is saying you should consider the "assignment" in the loop an equation. Call that equation eq[i,j], say.
Forget about doing a loop.
Make the set of equations S:={seq(seq( eq[i,j], i=1..N),j=0..N)};
Use solve to find the unknown u[i,j].
Give us the code as text or upload a worksheet if you need further assistance.

Since you could solve the problem exactly, dsolve ought to be able to do it too.
Following your lead dsolve can do it stepwise using that the first ode only has f[1], the second only f[1] and f[2], etc.
 

restart;
N:=4; alpha:=5*3.14/180; r:=10; Ha:=5; H:=1;
Rf:=diff(f[m-1](x),x,x,x)+2*alpha*r*sum(f[m-1-n](x)*diff(f[n](x),x),n=0..m-1)
+(4-Ha)*(alpha)^2*diff(f[m-1](x),x); 
ode:=diff(f[m](x),x,x,x)-CHI*diff(f[m-1](x),x,x,x)=h*H*Rf;
CHI:=`if`(m>1,1,0);
f[0]:=x->1-x^2;
## Using a list to keep the order:
sysL:=[seq(ode,m=1..N)];
bcsL:=[seq({f[m](0)=0,D(f[m])(0)=0,f[m](1)=0},m=1..N)];
for m from 1 to N do
  res[m]:=dsolve(eval({sysL[m]},{seq(res[i],i=1..m-1)}) union bcsL[m])
end do;
F:=f[0](x)+add(rhs(res[m]),m=1..N);
plot(eval(F,h=-0.9),x=0..1);
plot(eval(diff(F,x),x=1),h=-1.5..-0.2);

The plots look like the plots done with the numerical solution.

I tried your system and ended up using numerical solution.
 

restart;
N:=4; alpha:=5*3.14/180; r:=10; Ha:=5; H:=1;
Rf:=diff(f[m-1](x),x,x,x)+2*alpha*r*sum(f[m-1-n](x)*diff(f[n](x),x),n=0..m-1)
+(4-Ha)*(alpha)^2*diff(f[m-1](x),x); 
####
ode:=diff(f[m](x),x,x,x)-CHI*diff(f[m-1](x),x,x,x)=h*H*Rf;
CHI:=`if`(m>1,1,0);
f[0]:=x->1-x^2;
#### The ode system for f[i], i=1..N
sys:={seq(ode,m=1..N)};
bcs:=`union`(seq({f[m](0)=0,D(f[m])(0)=0,f[m](1)=0},m=1..N));
#### Introducing F as the sum avoiding using f (avoiding the possible confusion with f[i])
eq:=F(x)=add(f[m](x),m=0..N);
#### Since D(F)(1) is needed we introduce a second order ode for F:
odeF:=diff(eq,x,x);
#### Initial conditions for F:
icF0:=eval[recurse](eq,{x=0} union bcs);
eq1:=convert(diff(eq,x),D);
icF1:=eval[recurse](eq1,{x=0} union bcs);
#### The total system
SYS:=sys union {odeF}:
BCS:=bcs union {icF0,icF1}:
#### To obtain the two plots we let res return a list of two results:
res:=proc(h1) local sol,F1;
  if not h1::realcons then return 'procname(_passed)' end if;
  sol:=dsolve(eval(SYS,h=h1) union BCS,numeric,output=listprocedure);
  F1:=subs(sol,diff(F(x),x))(1);
  [sol,F1]
end proc;
#### Testing
res(-0.9);
#### The plots:
plots:-odeplot(res(-0.9)[1],[x,F(x)],0..1);
plot(res(h)[2],h=-1.5..-0.2,caption="D(F)(1) as a function of h");

The attempt to solve sys with bcs exactly fails:
 

dsolve(sys union bcs); #No answer

 

This happens if you use 2D input but doesn't if you use 1D input (a.k.a. Maple input).
I tried both in Maple 2018.1 and in Maple 2017.3.
No problem when using 1D input, but with 2D input I got the error message

Error, invalid loop statement termination

Yet another irritating difference between 1D input and 2D input.

MaplePrimes18-07-11_2D_problem.mw

Note. Just to make sure the 2D input problem has nothing to do with the fact that the statement doesn't do anything I tried this version:
 

if true then
  for i from 1 to 2 do
    i
  end do
else
  6
end if;

The same conclusion: Fine in 1D, error in 2D.

For numeric solution your system cannot contain vectors. You must rewrite the vector system so you have a system of scalar equations. That ought to be easy.
For help on this you must supply either code in text form or a worksheet (upload using the fat green arrow in the MaplePrimes editor). 

When using interface(typesetting=standard, the default in Maple 2015) the printing is done without parentheses.
In Maple 2018 the default setting of typesetting is extended.
So in Maple 2018 and Maple 2015 try:
 

restart;
interface(typesetting=extended); # The default in Maple 2018
t := table([1=table([a=123]) ]);
save t,"F:/MapleDiverse/test.m";
####
restart;
interface(typesetting=standard); # The default in Maple 2015
read "F:/MapleDiverse/test.m";
t[1][a] := 321;

In the first part (the saving) you can also try using interface(typesetting=standard) instead. It doesn't make any difference for the second part as long as that is done with typesetting=standard.

So my conclusion is that it just is a typesetting matter.

There has been a change from Maple 2015 to Maple 2016.
I tried the following in Maple 12, 17, 2015, 2016, 2017, and 2018:
 

restart;
type(fsolve,`module`);

In Maple 2016 and later fsolve is a module having ModuleApply as a local. Thus fsolve can be used as if it were a procedure.
In earlier version fsolve was just a procedure.
Incidentally, didn't you mean reinitialize instead of initialize?
If you try forget(fsolve,xxx=true) you don't get a complaint, thus a wrong optional argument of the type name= anything is just ignored.
The help page for forget has this statement about reinitialize:

"reinitialize = true or false (default: true)
This option only applies to the case when f is a procedure which is not a member of a module. In this case, if reinitialize=true, the corresponding procedure f is completely discarded and reloaded from its Maple archive, causing its remember table to be restored to its initial state." (my italics).

Try this:
 

op(2,eval(forget)); # forget is also a module
showstat(forget::ModuleApply);

 

Remove the range argument to dsolve and then do SOL(-2);
That keeps you from getting confused.
With the range argument just do SOL(-2);  before you go on.

As far as I can see your boundary value problem has no solution.
 

restart;
L:=75/10^6; # rational
eq1:=0.2079268293e-3*(diff(psi(x), x, x))-101.2059621*(diff(u(x), x, x))-101.2059621*(diff(w(x), x))*(diff(w(x), x, x));
##
eq2:=0.2079248162e-3+29.57317072*psi(x)-29.57317072*(diff(w(x), x))-1.996189024*10^(-9)*(diff(psi(x), x, x))-1.996189024*10^(-9)*(diff(w(x), x, x, x));
##
eq3:=151.8089432*(diff(w(x), x))^2*(diff(w(x), x, x))+101.2059621*(diff(u(x), x, x))*(diff(w(x), x))+101.2059621*(diff(u(x), x))*(diff(w(x), x, x))-0.2079268293e-3*(diff(psi(x), x, x))*(diff(w(x), x))-0.2079268293e-3*(diff(w(x), x, x))*(diff(psi(x), x))+29.57317070*(diff(w(x), x, x))-29.57317072*(diff(psi(x), x))-1.996189024*10^(-9)*(diff(psi(x), x, x, x))-1.996189024*10^(-9)*(diff(w(x), x, x, x, x))-p*(diff(w(x), x, x));
sys:=convert({eq1,eq2,eq3},rational): # for casesplit
bcs:=psi(0) = 0, psi(L) = 0, D(psi)(0) = 0, w(0) = 0, w(L) = 0, D(w)(0) = 0, D(w)(L) = 0, u(0) = 0, u(L) = 0;
with(DEtools):
res:=casesplit(sys,[u,w,psi]):
nops([res]);
for i from 2 to 5 do res[i] end do; #Special cases
#############################################################
## Testing the systems at x=0 and with bcs:
for i from 1 to 5 do
  temp:=op(1,res[i]);
  try
    R[i]:=eval[recurse](convert(temp,D),{x=0,bcs});
  catch "numeric exception":
    printf(StringTools:-FormatMessage(cat("i = %d\n",lastexception[2..-1],"\n")),i)
  end try;
end do:
R[2]; # 0 = 27224/3872072675, so no good
R[3]; # 0 = 27224/3872072675, so no good
R[4]; # 0 = 27224/3872072675, so no good
R[5]; # No contradiction
## Looking at the first two equations in res[1]:
temp:=op(1,res[1]):
R[1]:=eval[recurse](convert(temp[1..2],D),{x=0,bcs}); # 0 = 27224/3872072675, so no good
## So far the case res[5] is a possibility (and the only one)
res[5];
## Now look at that case 5 at x=L (and using bcs)
temp:=op(1,res[5]):
eval[recurse](convert(temp,D),{x=L,bcs}); # No contradiction
temp;
## We see that psi is constant (C), so odes for u and w are
sys5:=eval(temp[1..-2],psi(x)=C);
## We solve the w-equation:
bcsw:=op(select(has,{bcs},w));
## Notice that there are 4 boundary conditions
## One of those, w(L)=0, will be used to determine C.
solw:=dsolve({sys5[2],w(0)=0,D(w)(0)=0,D(w)(L)=0}); #Solving exactly
#Now determine C:
eval[recurse](solw,{x=L,w(L)=0});
C_result:=solve(%,{C}); # psi(x) = C
## This is already contradicting psi(L) = 0 which is required in bcs.

 

First 9 10 11 12 13 14 15 Last Page 11 of 150