## 11796 Reputation

15 years, 335 days

## OK in Maple 8...

I don't have Maple 7 on this machine, but I do have Maple 8.
I see all 4 lines in all 3 displays.
Since it appears that things are different in Maple 7 you could try using a list rather than a set. A list respects the order.
Even a sequence should work (it does in Maple 8).
It may be a matter of one plot covering the other?
So try permutations of:

`plots[display]( [lin[1],lin[2],lin[3],lin[4], squ,poinC,tetA], scaling=constrained,axes=normal);`

A slightly modified version of your code that works in Maple 2018.1 too:

```restart:
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Areas of portions of a square
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
with(plots):
with(plottools);
side:=6:A1:=16:A2:=20:A3:=32:
# Pick a point in the square: say (x1, y1)
poinC:=point([side/4,side/3]):

poinC:=point([side/4,side/3]);
cpoin:=convert(op(1,poinC),list);
# x1:=cpoin[1];
# y1:=cpoin[2];
lin[1]:=line([side/2,0],cpoin, color=black):
lin[2]:=line([0,side/2],cpoin, color=black):
lin[3]:=line([side/2,side],cpoin, color=black):
lin[4]:=line([side,side/2],cpoin, color=black):
#textA:=textplot([(cpoin[1]+side/2)/2,1,`A1`],align={ABOVE,RIGHT}):
#textA:=textplot([(1.5+side/2)/2,1,`A1`]):
tetA:=textplot([1,1.3,`A1`]):

squ:= polygon([[0,0], [0,side], [side,side],[side,0]], color=white):
# Plot without tetA
plots[display]({seq(lin[i],i=1..4),squ,poinC}, scaling=constrained,axes=normal);
# Plot with tetA
plots[display]({seq(lin[i],i=1..4),squ,tetA,poinC}, scaling=constrained,axes=normal);
# Plot with tetA (placed  differently)
plots[display]({squ,poinC,tetA,lin[1],lin[2],lin[3],lin[4]}, scaling=constrained,axes=normal);
```

## pdsolve/numeric...

pdsolve/numeric even insists that boundary conditions must be given on a region defined by lines parallel to the axes. Thus you cannot have a boundary condition like sh(H(t), t) = 0, unless H(t) is a constant, i.e. even if H(t) is actually known (e.g. H(t)=t), it wouldn't be allowed.

## An ode approach...

I cannot tell you what goes wrong. The logarithmic terms could possibly become imaginary at some point?
Anyway your pde1 can be considered an ode with the parameters f1(x), f2(x), f3(x), f4(x) as is done below.
The ugly and weird notation `&sigma;h3` and `&sigma;v3` is the 1D rendering of your 2D code. I hate 2D input with a passion. You see why.

```pde1;
ode:=subs(Lambda(x,t)=Lambda(t),f1(x)=F1,f2(x)=F2,f3(x)=F3,f4(x)=F4,pde1);
res:=dsolve({ode,Lambda(1.8936825887327868318*10^15) = F4},numeric,parameters=[F1,F2,F3,F4]);
## Testing. First setting parameters. We use x = 1000 here:
res(parameters=[F1=`&sigma;h3`(1000),F2=`&sigma;`(1000),F3=Jir3(1000),F4=Lambda3(1000)]);
res(Tf); #result
## A procedure to do that:
LL:=proc(xx,tt) if not type([xx,tt],list(realcons)) then return 'procname(_passed)' end if;
res(parameters=[F1=`&sigma;h3`(xx),F2=`&sigma;v3`(xx),F3=Jir3(xx),F4=Lambda3(xx)]);
subs(res(tt),Lambda(t))
end proc;
## Testing the procedure:
LL(1000,Tf);
## Plotting
plot(LL(1000,t),t=Ts..Tf); # x =1000
plot(LL(x,(Ts+Tf)/2),x=Xvp..Xp); # t= (Ts+Tf)/2
## The animation takes a while and is not very illuminating as it is:
plots:-animate(plot,[LL(x,t),x=Xvp..Xp],t=Ts..Tf);
```

Surely this couldn't be the most efficient way of producing a solution.

## Yes, an example...

Yes, I give you a very simple example where the initial and boundary conditions are defined by numerical procedures from dsolve/numeric:

```restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x,x);
ibcs:={u(x,0)=f(x),u(0,t)=g1(t),u(1,t)=g2(t)};
#Test with simple conditions
test_res:=pdsolve(pde,eval(ibcs,{f(x)=x,g1(t)=0,g2(t)=t+1}),numeric);
test_res:-animate(t=3);
# Now an odesys producing in fact the same functions f, g1, and g2:
odesys:={diff(f(x),x)=1,diff(g1(x),x)=0,diff(g2(x),x)=1,f(0)=0,g1(0)=0,g2(0)=1};
oderes:=dsolve(odesys,numeric,output=listprocedure);
F,G1,G2:=op(subs(oderes,[f(x),g1(x),g2(x)]));
res:=pdsolve(pde,eval(ibcs,{f=F,g1=G1,g2=G2}),numeric);
res:-animate(t=3);
test_res:-value()(.1234,3);
res:-value()(.1234,3);
```

## A bug...

Yes, this is a bug. The "solution" found by Maple 2018.1 (and also by Maple 2017.3) is
sol := u(t, x, y) = exp(-t*(x^2+2*x*y+y^2+1))
but it is wrong. It satisfies the initial condition, but not the pde.

In earlier versions (Maple 2016 and earlier) no solution was found.

I shall submit a bug report (aka an SCR).

## A bug...

This is clearly a bug in pdsolve/numeric.
A workaround for your case is to use the following pde (PDE2), which actually on the given range is identical to yours:

`PDE2:=PDE+(0=diff(L1(x,t),x)*Heaviside(-x));`

If you try

`simplify(PDE2) assuming x>0;`

you get PDE, but don't use that simplified version, keep PDE2 and use IBC := {L1(Xp, t) = 1, L1(x, Ts) = 1}.
Here is an extremely simple example of this same bug:

```restart;
pde:=diff(u(x,t),t)=u(x,t);
ibcs:={u(x,0)=x,u(0,t)=0};
pdsolve(pde,ibcs,numeric,time=t,range=0..1);
ibcs2:={u(x,0)=x};
pdsolve(pde,ibcs2,numeric,time=t,range=0..1);
sol:=pdsolve({pde,u(x,0)=x}); #Exact version
##### The trick:
pde2:=pde+(0=Heaviside(-x-1)*diff(u(x,t),x));
simplify(pde2) assuming x>=0;
ibcs:={u(x,0)=x,u(0,t)=0};
res:=pdsolve(pde2,ibcs,numeric,time=t,range=0..1);
res:-plot(t=5);
res:-animate(t=0..5);
## The same for the exact solution
plots:-animate(plot,[rhs(sol),x=0..1],t=0..5);
```

I think I may have sent an SCR (bug report) a long time ago about this, but I shall do it again just in case.

## A bug...

Clearly the answer from RealDomain:-solve({y=y,y<>0},y) is wrong and reveals a bug in RealDomain:-solve.

I think it can be traced from RealDomain:-solve via
SolveTools:-PolynomialSystemSolvers:-RealDomainFrontEnd
and
SolveTools:-PolynomialSystemSolvers:-RCRealSolver
to
SolveTools:-SemiAlgebraic,
which is being used here like this:
SolveTools:-SemiAlgebraic([0, y],{y});
Indeed

```SolveTools:-SemiAlgebraic([0, y<>0],{y});

```

returns [[y < 0], [0 < y]], which in the real domain is equivalent to y<>0.

## Solving for the highest derivatives...

@9009134 For your reformed program the system sys obtained by

`sys,bcs:=selectremove(has,dsys3,diff);`

still cannot immediately be solve for the highest derivatives hdiff:={diff(theta(x), x, x), diff(u(x), x, x, x, x), diff(w(x), x, x, x, x)}.
But sys[1] can be written
diff(u(x), x, x)+diff(w(x), x)*diff(w(x), x, x) = 0.
If that equation is replaced by the same differentiated twice the resulting new system can be solved for the highest deivatives hdiff.
Since you differentiated sys[1] twice you have 2 constraints:

diff(u(x), x, x)+diff(w(x), x)*diff(w(x), x, x) = 0 and
diff(u(x), x, x,x)+diff(w(x), x,x)^2 + diff(w(x), x)*diff(w(x), x, x,x) = 0.

Since the lhs of these will be constant you just need to evaluate them at one of the endpoints, e.g. 0
This gives you two new boundary conditions besides the 10 you already got. Thus you must remove two of the old ones.

Actually solving the new system for hdiff reveals that the denominator
(-1+400000*w(x))^6*(2*diff(u(x), x)+diff(w(x), x)^2)
appears.
This means that D(u)(0) and D(w)(0) cannot both be zero if you want to avoid a singularity at x=0. It's the same at x=0.001.
Thus in revising your bcs you should have this in mind.

## Solution is obtained in Maple 2016...

A solution is indeed found in Maple 2016.2. It takes a while and the result is so enormous that it is not printed, but a message is printed:
`[Length of output exceeds limit of 1000000]`

The result is there, but you should be happy that it is not printed.
Suppose you save the solution in sol.
Then try e.g.

```sol:=solve({e1, e2, e3, e4, e5}, {S__1, S__2, S__3, S__4, S__5});
nops([sol]); # 2
sol[1]; #Very simple
length(sol[2]); # 6216065
## Setting the parameters to 1 all of them for ease of testing:
params:=indets({e1, e2, e3, e4, e5},name) minus {S__1, S__2, S__3, S__4, S__5};
evalf[50](eval(sol[2],params=~1));
```

## The second is correct!...

Actually it is the second version that is correct. The first is false.

```restart;
eqns:=
[x=1, -1/(exp(h)+1/exp(h))^(1/2)/(exp(h)-1/exp(h))^(1/2)*(exp(h)^2-2+1/exp(h)^2)^(1/2)*x = tanh(h)^(1/2)];
eval(eqns,{x=1,h=2});
evalf(%); # Not equal
```

For h>0 eqns has no solutions for x. For h<=0 it does. (And nobody says that h couldn't be complex either).
Looking a little closer:

```eq2:=eval(eqns[2],x=1);
rhs(eq2)^2-lhs(eq2)^2;
convert(%,exp); # 0
```

It may be that solve in getting the first incorrect result squares eq2 as I just did.

## Example...

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
```

## Blasius' equation...

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
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.

## Difference solves homogeneous pde...

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

```pdetest(good,eqn_sys); # OK
##
HOM:=D[1, 1](y)(t, x)+y(t, x) = 0;
pdetest(solhom,HOM); # OK

```

## Two simple odes out of 3 for 2 unknowns...

Your question may be related to your question of June 19 2018:

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.

## Regression...

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.

 First 7 8 9 10 11 12 13 Last Page 9 of 149
﻿