Preben Alsholm

11796 Reputation

22 Badges

15 years, 335 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

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.

 

Maple 2018 has the ability to encrypt procedures. In "What is new"  it says:
"Maple 2018 lets you encrypt procedures. An encrypted procedure acts and behaves like any normal procedure in Maple, but it can not be viewed or debugged. This provides a way for professional users to share executable procedures while protecting their IP."

The result returned by int with fc is wrong.
 

restart;
f := ((1/2-I*t)^(-s)-(1/2+I*t)^(-s))/(2*I);
fc:=evalc(f);
res1:=int(f,t=0..infinity) assuming s>1;
int(fc,t=0..infinity) assuming s>1;
res2:=simplify(%);
# tests:
eval([res1,res2],s=2); # result [2,8]
evalf(Int(eval(fc,s=2),t=0..infinity)); # 2
evalf(Int(eval(f,s=2),t=0..infinity)); # 2
## Seeing all results:
int(fc,t=0..infinity,method=_RETURNVERBOSE) assuming s>1;
int(f,t=0..infinity,method=_RETURNVERBOSE) assuming s>1;

So for fc the "successful" methods are ftoc and ftocms and they agree and are wrong.
For f the "successful" methods are ftoc and ftocms as for fc, but in addition the correct meijerg.
## NOTE added: Today I cannot reproduce the wrong ftoc and ftocms results for f itself. They now  agree with the correct  meijerg result.

The simple exceptional solution not found by dsolve(ode) is found by `dsolve/IC/easy` when dsolve is given the initial value problem.
It can be found from the generic initial value problem by taking a limit:
 

restart;
ode:=diff(y(x),x)=x*ln(y(x)):
ic:=y(1)=y0: # generic ic
sol:=dsolve({ode,ic});
map(limit,sol,y0=1); # OK
## To see that the original initial value problem is handled by `dsolve/IC/easy`:
debug(`dsolve/IC/easy`);
dsolve({ode,y(0)=1});

 

Occasionally simplify makes an expression less simple than it was before or just not as simple as you had hoped.
In situations like that it can be worth while to try using map(simplify, expr).
Examples:
 

restart;
A := (1 - cos(s)^2)/sin(s)^2;
B := (1 - cos(t)^2)/sin(t)^2;
map(simplify,A*B);
##
S:=2*A+4*B;
simplify(S); # not simple
map(simplify,S); # 6

 

First 8 9 10 11 12 13 14 Last Page 10 of 149