## 11973 Reputation

16 years, 57 days

## Give the full path...

I made a .txt file and saved it at F:/test.txt.
Then I tried the example in the help page for readline, but replaced "infile.text" with the full path to my file.
Then no problem:

```restart;
file:="F:/test.txt";
while line <> 0 do
last := line;
end do;
```

## assuming odd and even...

Try this:

```f:=diff(ChebyshevT(n, r), r);
g:=simplify(eval(f, r = 0));
simplify(g) assuming n::odd;
simplify(eval(%,n=2*k+1)) assuming k::integer;
simplify(g) assuming n::even;
```

## Evolution equations...

@9009134 The numerical solver pdsolve/numeric can solve a system of pdes that in vector form can be written as
diff(u(x,t), t) = f(t, u(x,t), diff(u(x,t),x), .., diff(u(x,t),x\$n))
where u and f are vector valued.
If the given system is not first order in t it is converted to one (see the help page pdsolve/numeric).
Initial conditions of the form u(x,t0)=g(x) and boundary conditions must be given on two lines x=a and x=b.
The problems handled are thus time-based, i.e. the algorithm starts with t=t0 and u(x,t) = g(x) and then will work its way from there to any given t and along the way making sure that the boundary conditions are satisfied.

Your problem is not formulated that way. You have conditions given on the rectangle given by the lines x=0, x=L, z=0, z=h.
Your system could be rewritten (by pdsolve itself) as a system of the form given above if x is considered time, but then you cannot impose conditions on both of the lines x=0 and x=L. You must choose one of them, say x=0.  Furthermore since your system is second order in x (time) in both p1,p3, and Phi you must give initial conditions on the time derivatives too for x = 0.
A one dimensional analogue of the difference between your pde-bvp problem and the time based problem described here is the difference between a bvp problem for an ode and an initial value problem for the same ode.

### To make my comment clearer let me try doing what I said, i.e. consider x time and use z as the space variable. Notice my change of pd3. You probably meant one of the two occurrences of diff(Phi(x,z),x,x) to be diff(Phi(x,z),z,z). (If you meant what you wrote you must remove one of the boundary conditions for Phi).

```restart;
pd1 := -alpha*beta^2*diff(p1(x, z), x, x)-alpha*beta1^2*diff(p3(x, z), x, z)-alpha*beta2^2*diff(p1(x, z), z, z)-alpha*beta3*diff(p3(x, z), x, z)+alpha*p1(x, z)+diff(Phi(x, z), z) = f1(x)+z*f2(x);
pd2:=-alpha*beta1^2*diff(p1(x, z), x, z)-alpha*beta^2*diff(p3(x, z), z, z)-alpha*beta2^2*diff(p3(x, z), x, x)-alpha*beta3^2*diff(p1(x, z), x, z)+alpha*p3(x, z)+diff(Phi(x, z), z) = g(x);
##
## I have corrected your equation pd3:
## I changed one of the two occurrences of diff(Phi(x,z),x,x) to diff(Phi(x,z),z,z).
##
pd3 := -xi*(diff(Phi(x, z), x, x)+diff(Phi(x, z), z, z))+diff(p1(x, z), x)+diff(p3(x, z), z) = 0;
alpha := 1; beta1 := 2; beta2 := 3; beta3 := 4; h := 1; L := 5; xi := 6; beta := 4;
f1 := x->1/cos(x)^2;
f2 := x-> 1/(sin(x)^2+x^4);
g := x-> 1/cos(x)^4;
bcs:= Phi(x, 0) = 0,Phi(x, h) = 0, p1(x, 0) = 0, p1(x, h) = 0, p3(x, 0) = 0, p3(x, h) = 0;
ics:= Phi(0, z) = 0, p1(0, z) = 0, p3(0, z) = 0, D(Phi)(0,z)=0, D(p1)(0,z)=0, D(p3)(0,z)=0 ;
res:=pdsolve({pd1,pd2,pd3},{bcs, ics},numeric);
res:-plot3d(Phi,x=0..0.1);
```

Maybe it is helpful to show the simple conversion of your system to a system that it is of the form
diff(u(x,t), t) = f(t, u(x,t), diff(u(x,t),x), .., diff(u(x,t),x\$n)).

```restart;
pd1 := -alpha*beta^2*diff(p1(x, z), x, x)-alpha*beta1^2*diff(p3(x, z), x, z)-alpha*beta2^2*diff(p1(x, z), z, z)-alpha*beta3*diff(p3(x, z), x, z)+alpha*p1(x, z)+diff(Phi(x, z), z) = f1(x)+z*f2(x);
pd2:=-alpha*beta1^2*diff(p1(x, z), x, z)-alpha*beta^2*diff(p3(x, z), z, z)-alpha*beta2^2*diff(p3(x, z), x, x)-alpha*beta3^2*diff(p1(x, z), x, z)+alpha*p3(x, z)+diff(Phi(x, z), z) = g(x);
pd3 := -xi*(diff(Phi(x, z), x, x)+diff(Phi(x, z), z, z))+diff(p1(x, z), x)+diff(p3(x, z), z) = 0;
### Now do:
S1:={diff(p1(x,z),x)=q1(x,z),diff(p3(x,z),x)=q3(x,z),diff(Phi(x,z),x)=F(x,z)};
S2:=subs(S1,{pd1,pd2,pd3});
S2i:=solve(S2,{diff(q1(x,z),x),diff(q3(x,z),x),diff(F(x,z),x)});
### The system consisting of S1 and S2i together is now on the desired form:
map(print,S1 union S2i);
```

## Sets...

Some of your equations in eqs are sets. If you replace eqs by EQS as defined below you will get another error:

Error, (in dsolve/numeric/DAE/make_proc) number of unknown functions and equations must match, got 20 functions {C, C, C, C, C, C, Tg, Tg, Tg, Tg, Tg, Tg, Tg, Ts, Ts, Ts, Ts, Ts, Ts, Ts}, and 19 equations

```## Use this after your final definition of eqs:
EQS:=evalindets(eqs,set,op); # this removes the curly brackets
## Now use
sol := dsolve([op(EQS), op(ICs)], numeric, range = 0 .. tmax, abserr = 10^(-3), relerr = 10^(-3), stiff = true, maxfun = 0);
```

You must now think about the error message quoted.

The method as well as the notation used seem to be inspired by a Maple webinar by Sam Dao available from Maplesoft: Recorded webinars.

## A numerical solution...

Your odes are linear and the boundary value problems can probably all be solved exactly by dsolve in Maple directly. So these are obviously intended as test cases.
Kitonum solves the initial value problem u(0)=0, D(u)(0)=a exactly and use solve to find the exact value of a by requiring u(1)=0.
I shall give a numerical solution by a shooting method.
Here I shall use Problem 6:

```restart;
ode:=(1-x^2)*diff(u(x),x,x)-4*x*diff(u(x),x)-(1+x^2)*u(x)=x^2+cos(3*x);
sol:=dsolve({ode,u(-4)=8,u(-2)=-1}); #Exact solution directly
plot(rhs(sol),x=-4..-2);
## Now we find the numerical solution by shooting:
res:=dsolve({ode,u(-4)=8,D(u)(-4)=u1},numeric,parameters=[u1],output=listprocedure,abserr=1e-10,relerr=1e-9);
Up:=subs(res,u(x));
Q:=proc(u1) res(parameters=[u1]); Up(-2)+1 end proc;
plot(Q,0..10);
s:=fsolve(Q,0..20);
res(parameters=[s]); #Set the parameter u1 at s
plots:-odeplot(res,[x,u(x)-rhs(sol)],-4..-2); # The error
``` Finally it should be said that dsolve/numeric can solve the boundary value problem without using the shooting method. Just use

```res2:=dsolve({ode,u(-4)=8,u(-2)=-1},numeric);
```

######################################################
To get tabular results you can do as follows, where I use the exact solution sol and the approximate solution res (and concretely Up) found above.

```U:=unapply(rhs(sol),x); # The formula for the exact solution
N:=10: # A choice
interface(rtablesize=N+2);
X:=Vector([seq(-4+i*2/N,i=0..N)],datatype=float); # The x-values
Uexact:=evalf(U~(X)); # The exact solution at those x-values
Uapprox:=Up~(X); # The approximate solution  at those x-values
Udiff:=Uapprox-~Uexact; # The error: Approx-Exact
A:=<X|Uapprox|Uexact|Udiff>; # The data
L:=Array([x,u_approx(x),u_exact(x),u_approx-u_exact]); # Description
resM:=Matrix([[L],[A]]); # Final result``` ## 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,lin,lin,lin, 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;
# y1:=cpoin;
lin:=line([side/2,0],cpoin, color=black):
lin:=line([0,side/2],cpoin, color=black):
lin:=line([side/2,side],cpoin, color=black):
lin:=line([side,side/2],cpoin, color=black):
#textA:=textplot([(cpoin+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,lin,lin,lin}, 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 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 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; #Very simple
length(sol); # 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(eval(sol,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,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.

﻿