## 11796 Reputation

15 years, 335 days

## Singularity...

There could very likely be a singularity, not only a numerical problem.
Digits:=15:
# In the loop used this try..catch structure:

```try
sh0:=SH(0):
sv0:=SV(0):
Lx0:=LA(0):
print('i'=i);
catch:
print("Breaking",'i'=i); break
end try;
```

around the 3 assignments.
Furthermore, in the dsolve command I put in stiff=true, maxfun=10^6, and removed the range argument.
Then the loop ran just once. The break came at i = 2 probably because of an actual singularity strongly indicated here:

```sol(parameters);
odeplot(sol,[x,L(x)],4.238507..H-He);
odeplot(sol,[x,sh(x)],4.238507..H-He);
odeplot(sol,[x,sv(x)],4.238507..H-He);
odeplot(sol,[x,u(x)],4.238507..H-He);
```

## Execute again...

Looking at the output from the procedure definition of S1 (the blue stuff) I see psi inside the procedure body although the input appears to have phi (actually varphi).
My guess is that you had a first version defined like that. Then you realized that that was wrong and changed the input. You didn't execute the revised definition, however.

## Remove 'numeric'...

Just remove 'numeric':

```int(piecewise(0 < -k5*x1*x2+k1*x1 and -k5*x1*x2+k1*x1 < 1, 1, -k5*x1*x2+k1*x1 <= 0 and 1 <= -k5*x1*x2+k1*x1, 0)*piecewise(0 < -k5*x1*x2+k2*x2 and -k5*x1*x2+k2*x2 < 1, 1, -k5*x1*x2+k2*x2 <= 0 and 1 <= -k5*x1*x2+k2*x2, 0), k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1);
int(piecewise(0 < x+y and x+y < 1, 1, x+y <= 0 and 1 <= x+y, 0)*piecewise(0 < y*x^2 and y*x^2 < 1, 1, y*x^2 <= 0 and 1 <= y*x^2, 0), x = 0 .. 1, y = 0 .. 1);

That numeric integration takes a long time or has its problems is likely due to your integrands not being continuous.

## Nonsense...

Maybe it is a slight exaggeration to call your problem nonsense.
But my reason is this. Your system can in vector form be written as
X'(t) = f(t, X(t)), X(t0) = X0  (*)
where X and X0 are vectors, in your case of dimension 2, and f is a vector valued function.
The initial value problem (*) is equivalent to the integral equation
X(t) = X0 + int( f(s, X(s) ), s=t0 .. t);
As the value of an integral stays the same if the integrand is changed on a set of measure zero, your problem is equivalent to the simple problem

`sys:={diff(x(t),t)=y(t), diff(y(t),t)=-x(t), x(0)=0,y(0)=1};`

## A missing matching parenthesis...

In the definition of c a final parenthesis ')'  is missing. After that no complaint from the parser.

## Singularities...

Isn't it a problem that your function has infinitely many singularities on the positive real axis
(s[n], n=1..infinity) where s[n] -> infinity as n - > infinity
?
To see that there are singularities as described (and as expected because of the tan function) you could do:

```restart;
F:=1/(s^2*(1+(1-1/(2+2*s))/s)*(1+tan((1/2)*Pi/sqrt((1-1/(2+2*s))/s))/sqrt((1-1/(2+2*s))/s)));
F2:=normal(convert(F,sincos));
A:=op([1,1],indets(F2,specfunc(sin)))/Pi;
limit(A,s=infinity);
sp1,sp2:=solve(A=p,s);
evalf(seq(sp1,p=1..10));
dF2:=denom(F2);
seq(fsolve(dF2=0, s=eval(sp1,p=n)..eval(sp1,p=n+1)),n=0..9);
```

The sine term in the denominator dominates the cosine term when s > 1 and tends to zero as s->infinity. Thus the singularities will be close to those values of s that make the argument to sine an integral multiple of Pi. Those values are given by sp1 for p = 1..infinity.
The first 10 actual singularities are found in the last command by fsolve.

## Your example as an initial value problem...

Thanks to acer for suggesting major parts of the following approach.

Making use of dsolve 20000 times takes a while. So if the system could be turned into an initial value problem we can use the parameters option.

Although your system is really a boundary value problem since it has four functions with conditions given at H-He and one function (u) with u(0) = 0 we can consider instead posing the condition u(H-He)=0.
Your system is special. The ode for u has the form diff(u(x),x) = f( Jir(x) ) = g(x), i.e. it doesn't depend on u(x).
Thus the u(x) is just obtained by integration, u(x) = int(g(xi), xi = H-He .. x) (not that we shall do that literally).
The function u we actually want (ua, say) is then given by ua(x) = u(x)-u(0).
I have modified other parts, but that is really not important.

```restart;
Digits:=15:
phi0:=0.72:
rho0:=646.8:
E:=100*10^6:
nu:=1/3:
pc0:=0.3*10^6:
pv0:=0.6*10^6:
g:=9.80665:
Mp:=rho0*100/(10^6*365*24*60*60):
lambda:=E*nu/((1+nu)*(1-2*nu)):
mu:=E/(2*(1+nu)):
K:=lambda+2/3*mu:
eta:=100*10^6*10^6*365*24*60*60:
a:=1.15:
A:=-2/3+a*sqrt(3)/3:
B:=-1/3-a*sqrt(3)/3:
m:=5:
n:=3:
Ts:=60*10^6*365*24*60*60:
Tf:=2*Ts:
####
Md:=Mp*t:
phi:=1-(1-phi0)/Jir(x);
pc:=pc0*(ln(phi)/ln(phi0))^m;
mp:=-(m*(1-phi0)*pc)/(Jir(x)*phi*ln(phi));
Gp:=-(K+2*sqrt(3)/3*mu*a)/(K+a^2*mu+mp);
beta11:=K-2/3*mu;
beta21:=K+4/3*mu;
beta12:=(K-sqrt(3)/3*mu*a)*Gp+beta11;
beta22:=(K+2*sqrt(3)/3*mu*a)*Gp+beta21;
####
Delta1:=1/(rho0*g);
H:=beta21*Delta1*(1-exp(-g/beta21*Md));
Lambda:=1-rho0*g/beta21*(H-x);
####
sh:=beta11*ln(Lambda);
sv:=beta21*ln(Lambda);
Le:=evalf(exp(-pc0/(K+2*sqrt(3)/3*mu*a)));
Jire:=1;
He:=beta21*Delta1*(1-Le);
Te:=fsolve(H=He,t);
she:=eval(sh,{x=0,t=Te});
sve:=eval(sv,{x=0,t=Te});
ue:=-Mp*g/beta21*He;
unassign('Lambda','sh','sv','H'):
###############################
sys_ode:=diff(sv(x),x)=rho0*g/Lambda(x),diff(sh(x),x)=beta12/beta22*rho0*g/Lambda(x),diff(Lambda(x),x)=rho0*g/beta22,diff(u(x),x)=-Mp*g/beta22,diff(Jir(x),x)/Jir(x)=-Gp*diff(Lambda(x),x)/Lambda(x);
### Your ics has u(0)=0. The new ones:
icsNew:={sv(H-He)=sve,sh(H-He)=she,Lambda(H-He)=Le,Jir(H-He)=Jire,u(H-He)=0};
resNew:=dsolve({sys_ode} union icsNew,numeric,[sv(x),sh(x),Lambda(x),Jir(x),u(x)],parameters=[H],
output=listprocedure,abserr=1e-14,relerr=1e-12); #optimize may help a little
SV,SH,JIR,U:=op(subs(resNew,[sv(x),sh(x),Jir(x),u(x)]));
resNew(parameters=[He]); #Setting the parameter H to He for a test:
SV(0),SH(0),JIR(0),U(0);
#######
nn:=20000:
CC1:=Matrix(datatype=float):
CC2:=Matrix(datatype=float):
TR1:=Matrix(datatype=float):
TR2:=Matrix(datatype=float):
EN1:=Matrix(datatype=float):
EN2:=Matrix(datatype=float):
####
t0:=time():
dt:=(Ts-Te)/nn:
t:=Te:
HH:=He:
for i from 1 to nn+1 do
resNew(parameters=[HH]); #Setting the parameter H to HH
CC2(i,1..2):=Array([t/(10^6*365*24*60*60),HH]):
sh0:=SH(0):
sv0:=SV(0):
TR2(i,1..2):=Array([-1/3*sh0-2/3*sv0,evalf(sqrt(3)/3*(sh0-sv0))]):
Jir0:=JIR(0):
phix0:=1-(1-phi0)/Jir0:
EN2(i,1):=pc0*(ln(phix0)/ln(phi0))^m:
up:=U(HH-He)-U(0); # Since U is just an integral.
uh:=ue+up:
fvp:=evalf(A*sh0+B*sv0-(pv0*(ln(phix0)/ln(phi0))^n)):
if fvp>=0 then
break
else
dH:=dt*(Mp/rho0+uh):
HH:=HH+dH:
t:=t+dt:
end if:
end do:

time()-t0;
i;
```

The time used on my computer was 3.140s and i evaluates to 742.
MaplePrimes18-08-21_odesys_inits.mw

## _MonteCarlo...

Using method= _MonteCarlo and given a large enough epsilon:

`int(add(cat(x,i),i=1..10)^2,seq(cat(x,i)=0..1,i=1..10)],numeric,method=_MonteCarlo,epsilon=0.5e-4);`

I got 25.83247877.

## A dsolve/numeric workaround...

First an observation: In Maple 8 the 3 results of the evalf command in the following have no imaginary parts at all:

```M:=MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi);
seq(evalf(eval(M,x=xx)),xx=[0.4,0.5,0.6]);
```

In Maple 12, 15, 16, ..., 2018 there are very small imaginary parts in all 3.
This explains the difference in the results of this workaround for the integral:

```u:=(-5*ln(x)^4*Pi^4-20*ln(x)^2*Pi^4-8*Pi^4+120*MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi))/(120*Pi^4*(-1+x)^2);
ode:=diff(y(x),x)=u;
res:=dsolve({ode,y(4/10)=0},numeric); #OK in Maple 8, but not in 12, 15...2018
##A workaround for the workaround in 12, 15...2018
res:=dsolve({ode,y(4/10)=1e-22*I},numeric);
res(6/10);
```

The result of res(6/10) in Maple 2018.1 is -0.555168575558567e-2+1.00000000000001*10^(-22)*I
with the default settings of abserr and relerr.

## Bug?...

I tried the following and lost the kernel connection.

```restart;
u:=(-5*ln(x)^4*Pi^4-20*ln(x)^2*Pi^4-8*Pi^4+120*MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi))/(120*Pi^4*(-1+x)^2);
evalf(u);
int(%, x = 4/10 .. 6/10, numeric);
```

I shall report that as an SCR.

## Chop...

I once had that problem and then made a small procedure I called Chop lacking a better name.
Unlike evalf[n] which will turn exact expressions like Pi or sqrt(2) into floats with n digits Chop only works on existing floats.

```Chop:=proc (x::anything, n::posint) local op1, op2, k;
if not hastype(x, float) then return x end if;
if type(x, float) then
op1, op2 := op(x);
k := length(op1);
evalf(round(op1*10^(-k+n))*10^(op2-n+k), n)
elif type(x, complex(float)) then
procname(Re(x), n)+I*procname(Im(x), n)
else
evalindets(x,float,procname,n)
end if
end proc;

```

Examples:

```Chop(123.456789,4);
Chop(1.23456789,4);
evalf(f^Pi);
evalf(%); #Actually doesn't do anything!
Chop(%,4);
```

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

﻿