## 351 Reputation

13 years, 220 days

## Another (perhaps faster) approach -Simul...

The original problem has only 3 equations. I rewrote the model to write 3 equations for every value of delta (I used 21 data points for delta, you can increase it, but it might takDownload SimultaneousSolution.mwe more RAM).

Hence solving for 3*21 = 63 equations (21 different values of delta) seems to give faster results. Following CPU time is reported (includes time taken to dsolve and one plot. You should be able to return data till t= 0.5 and make 3D plots).

memory used=0.82GiB, alloc change=0 bytes, cpu time=5.84s, real time=5.85s, gc time=1.09s

If you change the model to real part and imaginary part, I expect dsolve to be much faster and RAM efficient.

PS., if anyone edits this, please avoid using ~/, shortcuts so that code is easy to follow and reverse compatible with older versions of Maple. (I am willing to learn shortcuts if gives me gain in cpu time/efficiency/RAM usage).

I am not able to display contents from my code. So, copying and pasting relevant lines of codes below.

restart;
Digits:=15;
15

phi:=0:lambda:=0.1:N:=5:M:=sqrt(N*(N+1))*exp(I*phi):omegap:=10:
dsys:=[diff(n(t),t)=-2*(n(t)-N)+(u(t)-M)*exp(-2*I*omegap*t/lambda)
+((z(t)-conjugate(M))*exp(2*I*omegap*t/lambda)),
diff(u(t),t)=-2*(1-I*delta)*u(t)
+2*(n(t)-N)*exp(2*I*omegap*t/lambda)+2*M,
diff((z(t),t))=-2*(1+I*delta)*z(t)
+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)
+2*conjugate(M)];
deltavals:=[seq(i-10.,i=0..20)];
[-10., -9., -8., -7., -6., -5., -4., -3., -2., -1., 0., 1., 2.,

3., 4., 5., 6., 7., 8., 9., 10.]

Ntot:=nops(deltavals);
21

vars:=[n(t)=NN[i1](t),u(t)=UU[i1](t),z(t)=ZZ[i1](t)];
[n(t) = NN[i1](t), u(t) = UU[i1](t), z(t) = ZZ[i1](t)]

for i from 1 to Ntot do eqs[i]:=evalf(subs(vars,i1=i,delta=deltavals[i],dsys));od:
ics:=seq(eval(NN[j](0))=0,j=1..Ntot),seq(eval(UU[j](0))=0,j=1..Ntot),seq(eval(ZZ[j](0))=0,j=1..Ntot):

Dsystot:=seq(op(eqs[i]),i=1..Ntot):

soln:=dsolve({Dsystot,ics},numeric,compile=true):
ss:=soln(0.5):

plot([seq([deltavals[i],subs(ss,NN[i](t))],i=1..Ntot)],style=point);

A simple procedure is written to find the cpu time/RAM usage for find the dsolve solution and then creating data for plots. Once sol(0.5) is found, solution at t= 0.5,
Simul:=proc()
local soln1,ss1,i;
soln1:=dsolve({Dsystot,ics},numeric,compile=true):
ss1:=soln1(0.5):
plot([seq([deltavals[i],subs(ss,NN[i](t))],i=1..Ntot)],style=point);
end proc;

with(CodeTools:-Profiling):
CodeTools:-Usage(Simul());
memory used=0.82GiB, alloc change=0 bytes, cpu time=5.84s, real time=5.85s, gc time=1.09s

## still slow...

If you are used to old mws format, then the best you can do today is to use worksheet mode with 1D input.

Note that current worksheet mode (from Maple 10 I believe) is Java based, and is very slow compared to classic worksheet model.

You can run current 32 bit versions in classic worksheet mode, but it crashes very frequently deleting all the commands. Dell seems to be better than hp computers for stability (in windows 10). Microsoft tablet is the worst (I could never install any stable mws format in that) and threw that laptop away.

## fsolve sucks big time...

If you have a good initial guess, try to write a simple Newton Raphson code.

As of now, fsolve will try to increase Digits arbitrarily and won't solve for nonlinear equations to a  specified tolerance. fsolve works only for nearly perfect residuals (= 0). For a larger system, it may not give any answers.

95 equations should be solvable. You can also try Optimization:-NLPSolve(1, []) with equations as constraints.

## Define a new and equivalent BVP such tha...

I don't recommend Shooting methods. In particular, single shooting method will fail for many nonlinear and stiff problems. Multiple shooting is better. Use of optimization (combining IVP or BVP solver with an optimizer) is inefficient and should be used only for optimal control or parameter estimation unless you know that the BVPs have a reasonable initial guess for unknown boundary conditions at one end.

Divide x = 0.. 2

as x = 0..1 (region 1) , and x2 = x-1 = 0..1 (region 2)

Define y as y in region 1 and y2 in region 2.

Without loss of generality replace x2 with x. Make sure you don't have any explicit functions of x (non-homogeneous non-constant terms in your equation. If there is any, use a dummy variable z and use dz/dx = 1)

We need additional boundary conditions. Derivatives and second derivatives are continuous at region 1/region 2 interface. So

dy/dx(1) = dy2/dx(0), etc.

For details on related approach (to scale any domain to 0 to 1), see past papers on BVP solvers or one of our past papers, (figure 1 in particular)

http://depts.washington.edu/maple/pubs/40_JES_OC_reformulation.pdf

You might want to see a different approach to solve BVPs (in particular if Maple's BVP solver fails or if you want to solve a largescale BVP).

https://www.mapleprimes.com/posts/208499-A-New-And-Efficient-Code-In-Maple-For

 > restart;
 > Digits:=15;
 (1)
 > sola:=dsolve([diff(y(x), x\$3)+diff(y(x), x\$2)+y(x)=1, y(0)=0, y(1)=0, y(2)=1], [y(x)]):
 > analoutput:=[evalf(subs(sola,x=0.5,y(x))),evalf(subs(sola,x=1.5,y(x)))];
 (2)
 > eq1:=diff(y(x), x\$3)+diff(y(x), x\$2)+y(x)=1;
 (3)
 > eq2:=subs(y(x)=y2(x),eq1);
 (4)
 > solnumeric:=dsolve([eq1,eq2,y(0)=0,y(1)=0,y2(0)=0,y2(1)=1,D(y)(1)=D(y2)(0),(D@@2)(y)(1)=(D@@2)(y2)(0)],[y(x),y2(x)],numeric):
 > s1:=solnumeric(0.5);
 (5)
 > numoutput:=[subs(s1,y(x)),subs(s1,y2(x))];
 (6)
 >

Dr. Venkat Subramanian

## Try differentiating, there are no DAE so...

Despite what Maple claims and what the help files say, there are no DAE solvers in Maple. Maple converts DAEs to dy/dt = f format and solves.

See my approach. See my past posts that helps avoid using fsolve. You have to make sure this code is correct( there may be more than one solution for initial condition).

In the future if possible please post y1, y2, etc instead of __t which took time for me to figure out.

 > restart;
 > Digits:=15;
 (1)
 > sys:={Q(0) = 0, Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2), diff(Q(t), t) = (0.1375e-1*47.6035070726347)*(T__1(t)-T__1s(t))*(T__1(t)+T__1s(t))^.438263122318020*((T__1(t)-T__1s(t))^.327228811371115), diff(Q(t), t) = (0.1375e-1*47.6035072491656)*(T__2s(t)-T__2(t))*(T__2(t)+T__2s(t))^.438263121701134*((T__2s(t)-T__2(t))^.327228811154163)};
 (2)
 > indets(sys);
 (3)
 > eqs:=[Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2), diff(Q(t), t) = (0.1375e-1*47.6035070726347)*(T__1(t)-T__1s(t))*(T__1(t)+T__1s(t))^.438263122318020*((T__1(t)-T__1s(t))^.327228811371115), diff(Q(t), t) = (0.1375e-1*47.6035072491656)*(T__2s(t)-T__2(t))*(T__2(t)+T__2s(t))^.438263121701134*((T__2s(t)-T__2(t))^.327228811154163)];
 (4)
 > eqs[4]:=subs(eqs[3],eqs[4]);
 (5)
 > eqs[5]:=subs(eqs[3],eqs[5]);
 (6)
 > eqae:=[eqs[1],eqs[2],eqs[4],eqs[5]];
 (7)
 > eqs1:=eval(subs(Q(t)=0,eqae));
 (8)
 > eqs2:=subs(T__1(t)=T1,T__2(t)=T2,T__1s(t)=T1s,T__2s(t)=T2s,eqs1);
 (9)
 > sol1:=fsolve({op(eqs2)});
 (10)
 > ics:=[Q(0)=0,T__1(t)=T1,T__2(t)=T2,T__1s(t)=T1s,T__2s(t)=T2s];
 (11)
 > ics1:=subs(sol1,t=0,ics);
 (12)
 > sys2:=diff(eqs[1],t),diff(eqs[2],t),diff(eqs[4],t),diff(eqs[5],t),eqs[3];
 (13)
 >
 > sol2:=dsolve({sys2,op(ics1)},type=numeric,implicit=true):
 > sol2(1);
 (14)
 > plots:-odeplot(sol2,[t,Q(t)],0..1);
 > plots:-odeplot(sol2,[t,T__1(t)],0..1);
 >

 > restart;
 > Digits:=15;
 (1)
 > sys:={Q(0) = 0, Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2), diff(Q(t), t) = (0.1375e-1*47.6035070726347)*(T__1(t)-T__1s(t))*(T__1(t)+T__1s(t))^.438263122318020*((T__1(t)-T__1s(t))^.327228811371115), diff(Q(t), t) = (0.1375e-1*47.6035072491656)*(T__2s(t)-T__2(t))*(T__2(t)+T__2s(t))^.438263121701134*((T__2s(t)-T__2(t))^.327228811154163)};
 (2)
 > indets(sys);
 (3)
 > eqs:=[Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2), diff(Q(t), t) = (0.1375e-1*47.6035070726347)*(T__1(t)-T__1s(t))*(T__1(t)+T__1s(t))^.438263122318020*((T__1(t)-T__1s(t))^.327228811371115), diff(Q(t), t) = (0.1375e-1*47.6035072491656)*(T__2s(t)-T__2(t))*(T__2(t)+T__2s(t))^.438263121701134*((T__2s(t)-T__2(t))^.327228811154163)];
 (4)
 > eqs[4]:=subs(eqs[3],eqs[4]);
 (5)
 > eqs[5]:=subs(eqs[3],eqs[5]);
 (6)
 > eqae:=[eqs[1],eqs[2],eqs[4],eqs[5]];
 (7)
 > eqs1:=eval(subs(Q(t)=0,eqae));
 (8)
 > eqs2:=subs(T__1(t)=T1,T__2(t)=T2,T__1s(t)=T1s,T__2s(t)=T2s,eqs1);
 (9)
 > sol1:=fsolve({op(eqs2)});
 (10)
 > ics:=[Q(0)=0,T__1(t)=T1,T__2(t)=T2,T__1s(t)=T1s,T__2s(t)=T2s];
 (11)
 > ics1:=subs(sol1,t=0,ics);
 (12)
 > sys2:=diff(eqs[1],t),diff(eqs[2],t),diff(eqs[4],t),diff(eqs[5],t),eqs[3];
 (13)
 >
 > sol2:=dsolve({sys2,op(ics1)},type=numeric,implicit=true):
 > sol2(1);
 (14)
 > plots:-odeplot(sol2,[t,Q(t)],0..1);
 > plots:-odeplot(sol2,[t,T__1(t)],0..1);
 >

## See code attached...

p = df/deta, so, p, theta and F are solved.

Finite difference discretization is used.

Second order (central difference) for first derivative gives more oscillations. So first order approximation is used.

Code not very stable as fsolve does not work. Even NLPSolve obj = 1 is not robust.

if you think these plots are meaningful. (sometimes P turns out to be negative), then you can build on this code for large values of inf (increae until there is no change) and increase N.

Check for typos and mistakes in the code as well.

 > restart;
 > Digits:=15;
 (1)
 > equ1 := (1+2*n)*f(eta)*(diff(theta(eta), eta))/(1+3*n)-(diff(theta(eta), `\$`(eta, 2))) = 0;
 (2)
 > equ2 := ((1+n)*(diff(f(eta), eta))^2/(1+3*n)-(1+2*n)*f(eta)*(diff(f(eta), eta, eta))/(1+3*n))/Bo+(diff(f(eta), `\$`(eta, 3)))^n-theta(eta) = 0;
 (3)
 > Bo:=1;n:=2;
 (4)
 > inf:=4;
 (5)
 > bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(inf) = 0, theta(0) = 1, theta(inf) = 0;#(D@@3)(f)(0) = z(0);
 (6)
 > equ1;equ2;
 (7)
 > equ2:=subs(diff(f(eta),eta)=p(eta),equ2);
 (8)
 > equ3:=diff(f(eta),eta)-p(eta)=0;
 (9)
 > N:=20;
 (10)
 > Eq1[0]:=T[0]=1;Eq2[0]:=P[0]=0;Eq3[0]:=F[0]=0;
 (11)
 > #for i from 1 to N do Eq1[i]:=subs(diff(theta(eta),eta\$2)=(T[i+1]-2*T[i]+T[i-1])/h^2,diff(theta(eta),eta)=(T[i+1]-T[i-1])/2/h,f(eta)=F[i],equ1);od;
 > for i from 1 to N do Eq1[i]:=subs(diff(theta(eta),eta\$2)=(T[i+1]-2*T[i]+T[i-1])/h^2,diff(theta(eta),eta)=(T[i]-T[i-1])/1/h,f(eta)=F[i],equ1);od;
 (12)
 > for i from 1 to N do #Eq2[i]:=subs(diff(p(eta),eta\$2)=(P[i+1]-2*P[i]+P[i-1])/h^2,diff(p(eta),eta)=(P[i+1]-P[i-1])/2/h,f(eta)=F[i],theta(eta)=T[i],p(eta)=P[i],equ2);od; Eq2[i]:=subs(diff(p(eta),eta\$2)=(P[i+1]-2*P[i]+P[i-1])/h^2,diff(p(eta),eta)=(P[i]-P[i-1])/1/h,f(eta)=F[i],theta(eta)=T[i],p(eta)=P[i],equ2);od;
 (13)
 > #for i from 1 to N do Eq3[i]:=subs(diff(f(eta),eta)=(F[i+1]-F[i-1])/2/h,p(eta)=P[i],equ3);od;
 > for i from 1 to N do Eq3[i]:=subs(diff(f(eta),eta)=(F[i]-F[i-1])/1/h,p(eta)=P[i],equ3);od;
 (14)
 > Eq1[N+1]:=T[N+1]=0;Eq2[N+1]:=P[N+1]=0; #Eq3[N+1]:=3*F[N+1]-4*F[N]+F[N-1]=0; Eq3[N+1]:=F[N+1]-F[N]=0;
 (15)
 > eqs:=seq(Eq1[i],i=0..N+1),seq(Eq2[i],i=0..N+1),seq(Eq3[i],i=0..N+1),h=inf/(N+1):
 > guess:=[seq(T[i]=0.9,i=0..N+1),seq(P[i]=-0.1,i=0..N+1),seq(F[i]=0,i=0..N+1),h=inf/(N+1)];
 (16)
 > #sol:=fsolve({eqs},{op(guess)});#fsolve may not work
 >
 > sol1:=Optimization:-NLPSolve(1,[eqs],initialpoint=guess);sol:=sol1[2];
 (17)
 > plot(subs(sol,[seq([i*h,T[i]],i=0..N+1)]),axes=boxed);
 > plot(subs(sol,[seq([i*h,F[i]],i=0..N+1)]),axes=boxed);
 > plot(subs(sol,[seq([i*h,P[i]],i=0..N+1)]),axes=boxed);
 >
 >
 >
 >
 >
 >
 >
 >

## FWIW...

My experience with optmization in Maple for the last 10 years or so.

Help pages are poorly documented. So, you have to try many things to figure out what you want to do.

NLP problems can be done with different approaches in Maple.

As of now, more reliable answers are arrived for sqp approach (bounded and constrained problems) only if you state the problem in algebraic form (i.e, variables y1 to yn, bounds, etc in explicit form). Depsite what the help page claims, the matrix vector f(Y) form or using functional form f(u1,u2..) is not reliable. There seems to be some  bug in finding the numerical jacobian or hessian that will make Maple suggest that "initial guess meets first order optimality conditions, etc). method = nonlinearsimplex is reasonably robust and can be used only for unbounded and uncostrained problems.

Also, GAMS has access to many sophisticated optimizers. Maple's nonlinear constrained optimization is limited to regular sqp. So, you will be limtied by the size of the problem you can handle. This is because, Maple stores the entire Jacobian, and does not use sparse storage or solver approaches. It is easy to verify. Just take d^y/dx^2 = y^2 with bcs, x =0, y =1 and x = 1, y =0. Convert to finite difference form, and give 1 as objective and give bounds of 0..1 to for y. You will hit the size limit N  (number of node points) very quickly (perhaps N = 200 or 500) is the maximum size of optimization variables Maple can handle.

I have also used the Global optimization package (there is no gurantee for global optimization for any commerical package, despite the name) currently offered. The current package includes Genetic algorithm, but the previous global optimization offered by Pinto (during Maple 12,13) was more robust and had a good reduced gradient approach which was fast.

Mathematica offers access to KNITRO which is very robust. I think GAMS provides access to SNOPT and KNITRO, may be even IPOPT. Of course, people routinely do MEX files from MATLAB for Ipopt, knitro, etc.

I wish Maple provides access to KNITRO, SNOPT and IPOPT. I will be surprised if it happens considering the direction Maple has taken recently (focussing more on elementary stuff, compared to high-end users). That is a shame as Maple has the ability to do analytic Jacobian.

In our group, we run IPOPT from Maple by writing the C files for objective, gradients, constraints, jacobians and hessians. Proper documentation and more useful help files for externalcalling will help me and many to create a dll or package that can be added to Maple in the future.

## check past post and explanation given by...

http://www.mapleprimes.com/questions/125273-Dsolve-Events-How-To-Control-For-A-Sign-Change#comment125661

## I have a solution...

Bendesarts

See below and attached.

The model you have might be from a Hamiltonian. This means that all the solvers in Maple will fail. Read about geometric numerical integration. Basically we need algorithms that conserve energy. If you know which one is like displacement, and which one is like velocity, you may be able to come up with symplectic Euler like explicit schemes.

Fortunately, Gauss implicit Runge Kutta works for these type of problems (symplectic integrators). The second order Gauss method is the implicit mid point method.

That is for a system of ODEs dY/dt = F, Y1 = Y0+h*F(Ymid) where Ymid is given by Ymid = (Y0+Y1)/2. If F is nonlinear, then one needs to do a Newton Raphson iteration at every time step. This is a A stable (not L stable ) method good for some stiff problems, but does not have stiff decay property. Also not algebraically stable.

More about the property of this method later. The code below works. You can see that both dt = 0.1 and dt = 0.2 give same results from 0 to 400 for time.

Few ways to improve the code. Let me know if there is interest

1. Write the code in hardware floats (Speed up). This will mean that fsolve should be avoided and Newton Raphson should be done without using LinearSolve or any linear algebra package. Read my past code/post on Newton Raphson. Convert the same to compiled form for even more speed. If not for speed, do this for robustness as fsolve will fail and crash for large systems of equations.

2. Write adaptive time stepping instead of constant dt. This wil speed up the process and will also guarantee that results meet the tolerance specified.

3. Use higher order methods  - Higher order algorithms or simply using extrapolation methods. This might compromise on stability (in particular extrapolation).

4. Write schemes other than implicit mid point to march forward in time.

5. Write this and related (Adatptive time stepping with events, hardware floats/compiled form) for generic systems of ODEs and DAEs. The code is posted for ease of learning, not optimized for efficiency or performance.

Hope this is useful

 > restart;
 > r:=sqrt((x(t)/a)^2+(z(t)/b)^2); eqx:=diff(x(t),t)=alpha*(1-r^2)*x(t)+w*a/b*z(t); eqz:=diff(z(t),t)=beta*(1-r^2)*z(t)-w*b/a*x(t); EqSys:=[eqx,eqz]; params := alpha=1, beta=1, a=0.4, b=0.2, w=1; EqSys := eval([eqx,eqz], [params]); xmax := 0.8; zmax := 0.4; tmax := 400; ic:=[x(0)=0.4, z(0)=0];
 (1)
 > DEplot(EqSys, [x(t),z(t)], t= 0..tmax, [ic],linecolor=black, thickness=1,x(t)=-xmax..xmax, z(t)=-zmax..zmax, scaling=constrained,arrows=none);
 (2)
 > sol:=dsolve({op(EqSys),op(ic)},type=numeric,range=0..400):
 > with(plots):
 > odeplot(sol,[x(t),z(t)],t=0..400,thickness=3,axes=boxed);
 > sol(400);
 (3)
 > r:='r':
 > r:=sqrt((x(t)/a)^2+(z(t)/b)^2); eqx:=diff(x(t),t)=alpha*(1-r^2)*x(t)+w*a/b*z(t);
 (4)
 > eqz:=diff(z(t),t)=beta*(1-r(t)^2)*z(t)-w*b/a*x(t);
 (5)
 > ic2:=x(0)=0.4,z(0)=0;
 (6)
 (7)
 > Digits:=15;
 (8)
 > eq1:=(lhs-rhs)(subs(diff(x(t),t)=(xnew-xold)/dt,diff(z(t),t)=(znew-zold)/dt,x(t)=(xnew+xold)/2,z(t)=(znew+zold)/2,eqx));
 (9)
 > eq2:=(lhs-rhs)(subs(diff(x(t),t)=(xnew-xold)/dt,diff(z(t),t)=(znew-zold)/dt,x(t)=(xnew+xold)/2,z(t)=(znew+zold)/2,eqz));
 (10)
 > ss:=proc(Y0,dt) local xold,zold,eq1,eq2,s2; xold:=Y0[1];zold:=Y0[2]; eq1:=(xnew-xold)/dt-(1-6.250000000*(1/2*xnew+1/2*xold)^2-25.00000000*(1/2*znew+1/2*zold)^2)*(1/2*xnew+1/2*xold)-1.00000000000000*znew-1.00000000000000*zold; eq2:=(znew-zold)/dt-(1-6.250000000*(1/2*xnew+1/2*xold)^2-25.00000000*(1/2*znew+1/2*zold)^2)*(1/2*znew+1/2*zold)+.250000000000000*xnew+.250000000000000*xold; s2:=fsolve({eq1,eq2},{xnew=xold,znew=zold}); Vector(2,[subs(s2,xnew),subs(s2,znew)]); end proc;
 (11)
 > Y[0]:=Vector(2,[0.4,0]);
 (12)
 > ss(Y[0],0.1);
 (13)
 > tmax:=400.;N:=4000;dt:=tmax/N;
 (14)
 > for i from 1 to N do Y[i]:=ss(Y[i-1],dt);od:
 > p1:=plot([seq([Y[i][1],Y[i][2]],i=0..N)],axes=boxed):
 > N:=2000;dt:=tmax/N;
 (15)
 >
 > for i from 1 to N do Y[i]:=ss(Y[i-1],dt);od:
 > p2:=plot([seq([Y[i][1],Y[i][2]],i=0..N)],axes=boxed,color=blue,thickness=3,axes=boxed):
 > display({p1,p2});
 >

## False transient and robust transient met...

Al86

See code attached below. y and z are used instead of x and y. phi is used instead of u. Basically artificial time derivatives are added and finite difference or other methods are used in x and y  in transient method of lines. For your problem, that may not work (if you have sqrt or stiff problems).

Attached code does both standard false transient and new a robust false transient method of lines we proposed. You can read more about it at http://depts.washington.edu/maple/falsetransient_new.html

Hope you find the approach useful. For higher number of node points, using stiff =true and compile =true will help. Stady state is reached in t = 0.05. I think the profile makes sense. You might want to increase the number of node points.

 > restart;
 > with(plots);
 (1)

Height of the domain (y-coordinate)

 > h:=1;
 (2)

Length of the domain (z-coordinate)

 > L:=1;
 (3)
 >

Governing Equation (Laplace's Equation)

 > Eq1:=diff(phi(y,z),y\$2)+diff(phi(y,z),z\$2)-sqrt(phi(y,z))-diff(phi(y,z),y)^2/phi(y,z)^(3/2);
 (4)

Boundary condition at y=0

 > bcy1:=phi(y,z)-1;#diff(phi(y,z),y);
 (5)

Boundary condition at y=1

 > bcy2:=diff(phi(y,z),y);
 (6)

Boundary condition at z=0

 > bcz1:=phi(y,z)-1;
 (7)

Boundary condition at z=1

 > bcz2:=diff(phi(y,z),z);
 (8)

Finite difference scheme to be used (2nd order-central difference for 2nd derivative)

 > d2phidy2:=(phi[m+1,n]-2*phi[m,n]+phi[m-1,n])/dely^2; d2phidz2:=(phi[m,n+1]-2*phi[m,n]+phi[m,n-1])/delz^2;
 (9)

Finite difference scheme to be used (2nd order-central difference for 1st derivative)

 > dphidy:=(phi[m+1,n]-phi[m-1,n])/(2*dely);
 > dphidz:=(phi[m,n+1]-phi[m,n-1])/(2*delz);
 (10)

3 Point forward and backward differences for 1st derivative (to be applied at boundary conditions)

 > dphidyf:=(-phi[2,n]+4*phi[1,n]-3*phi[0,n])/(2*dely); dphidyb:=(phi[Numy-1,n]-4*phi[Numy,n]+3*phi[Numy+1,n])/(2*delz); dphidzf:=(-phi[m,2]+4*phi[m,1]-3*phi[m,0])/(2*delz); dphidzb:=(phi[m,Numz-1]-4*phi[m,Numz]+3*phi[m,Numz+1])/(2*delz);
 (11)

Number of interior node points in y

 > Numy:=5;
 (12)

Number of interior node points in z

 > Numz:=5;
 (13)
 >
 > dely:=(h)/(Numy+1);
 (14)
 > delz:=(L)/(Numz+1);
 (15)
 >

Develop finite difference equations for z boundary conditions for all y

 > for j from 0 to Numy+1 do Eq[j,0]:=subs(diff(phi(y,z),z)=dphidzf,phi(y,z)=phi[j,0],m=j,bcz1): Eq[j,Numz+1]:=subs(diff(phi(y,z),z)=dphidzb,phi(y,z)=phi[j,Numz+1],m=j,bcz2): od;
 (16)

Develop finite difference equations for y boundary conditions for all z

 > for j from 0 to Numz+1 do Eq[0,j]:=subs(diff(phi(y,z),y)=dphidyf,phi(y,z)=phi[0,n],n=j,bcy1): Eq[Numy+1,j]:=subs(diff(phi(y,z),y)=dphidyb,phi(y,z)=phi[Numy+1,n],n=j,bcy2): od;
 (17)
 > #printlevel:=2;
 >
 >

Develop finite difference equations for all interior points using the governing equation

 > for j from 1 to Numy do for k1 from 1 to Numz do Eq[j,k1]:=subs(diff(phi(y,z),y\$2)=d2phidy2,diff(phi(y,z),z\$2)=d2phidz2,diff(phi(y,z),y)=dphidy,diff(phi(y,z),z)=dphidz,phi(y,z)=phi[j,k1],n=k1,m=j,Eq1); od;od;
 > #printlevel:=1;

Collect all equations into a single list

 > eqs:=[seq(seq(Eq[p,q],p=0..Numy+1),q=0..Numz+1)]:

Collect all variables into a single list

 > vars:=[seq(seq(phi[i,j],i=0..Numz+1),j=0..Numy+1)]:
 >

Count number of equations

 > n1:=nops(eqs);
 (18)
 >

Convert all variables from the form phi[i,j] to YY[i]

 > vars2:=[seq(vars[i]=YY[i](t),i=1..n1)];
 (19)

Perturbation parameter to be used

 > mu:=1e-3:

Substitute new variables into the equations

 > Eqs:=subs(vars2,eqs):

Eqs2 is the standard false transient formulation

 > Eqs2:=seq(diff(YY[i](t),t)=Eqs[i],i=1..n1):

Eqs3 is the perturbation approach described in the paper

 > Eqs3:=seq(mu*diff(Eqs[i],t)=-Eqs[i],i=1..n1):
 >

This is an initial guess for all values of phi to be used in the IVP solver

 > ics2:=seq(YY[i](0)=1,i=1..n1):

Solver the standard false transient formulation with Maple's dsolve

 > temp:=time[real]():sol2a:=dsolve({Eqs2,ics2},type=numeric,implicit=true,stiff=true,compile=true):time[real]()-temp;
 (20)

Solver the perturbation formulation with Maple's dsolve

 > temp:=time[real]():sol3a:=dsolve({Eqs3,ics2},type=numeric,implicit=true,stiff=true,compile=true):time[real]()-temp;
 (21)
 >
 >

Plot the convergence of the standard false transient (red) and the perturbation approach (green). YY[1] is phi at y=0, z=0

 > t11:=time[real]():p2:=odeplot(sol2a,[t,YY[25](t)],0..10):t11-time[real](); t11:=time[real]():p3:=odeplot(sol3a,[t,YY[25](t)],0..10,color=green):t11-time[real](); display(p2,p3);

 Warning, cannot evaluate the solution further right of 1.1718747, probably a singularity

Calculate converged value from false transient approach (value at pseudo time=50)

 > s2:=sol2a(50);
 Error, (in sol2a) cannot evaluate the solution further right of 1.1718747, probably a singularity

Calculate converged value from perturbation approach (value at pseudo time=10)

 > s3:=sol3a(10);
 (22)
 > odeplot(sol3a,[t,YY[9](t)],0..0.05);
 >

Convert the solution (from s3, as a list), as an array

 > for p from 0 to Numy+1 do
 > for q from 0 to Numz+1 do phi3[p,q]:=subs(vars2,s3,phi[p,q]); yt:=p/(Numy+1); zt:=q/(Numz+1); od; od:
 >
 >
 >

Plot the solution as found from the perturbation approach

 > sol3:=matrix([seq([seq(phi3[p,q],p=0..Numy+1)],q=0..Numz+1)]):
 > matrixplot(sol3,axes=boxed);
 >
 >
 >
 >
 >
 >
 >

## Try DAE form - there are 3 solutions, al...

See below. The problem is converted to index 1 DAE and solved. This problem might have multiple solutions (bifurcation) and might become index 2 DAE. I will explore more. Even at t= 0, there are 3 possible solutions.

 > restart:
 > odes:=diff(x1(t),t)*diff(x2(t),t\$2)*sin(x1(t)*x2(t))+5*diff(x1(t),t\$2)*diff(x2(t),t)*cos(x1(t)^2)+t^2*x1(t)*x2(t)^2=exp(-x2(t)^2),diff(x1(t),t\$2)*x2(t)+diff(x2(t),t\$2)*diff(x1(t),t)*sin(x1(t)^2)+cos(diff(x2(t),t\$2)*x2(t))=sin(t); ics:=x1(0)=1,D(x1)(0)=1,x2(0)=2,D(x2)(0)=2; odes2:=subs(diff(x1(t),t\$2)=yp2,diff(x2(t),t\$2)=yp4,diff(x1(t),t)=Y[2],diff(x2(t),t)=Y[4],x1(t)=Y[1],x2(t)=Y[3],{odes});
 (1)
 > odes3:=subs(diff(x1(t),t\$2)=z1(t),diff(x2(t),t\$2)=z2(t),diff(x1(t),t)=x1t(t),diff(x2(t),t)=x2t(t),{odes});
 (2)
 > odes4:=odes3 union {diff(x1(t),t)=x1t(t),diff(x1t(t),t)=z1(t),diff(x2(t),t)=x2t(t),diff(x2t(t),t)=z2(t)};
 (3)
 > Digits:=15;
 (4)
 > eqic:=eval(subs(x2(t)=2,x1(t)=1,x1t(t)=1,x2t(t)=2,z1(t)=z10,z2(t)=z20,t=0.,odes3));
 (5)
 > z10:=solve(eqic[1],z10);
 (6)
 > eqic[2];
 (7)
 > zs:=[solve(eqic[2],z20)];
 (8)
 > z10:=subs(z20=alpha,z10);
 (9)
 > sol:=dsolve({op(odes4),x2(0)=2,x1(0)=1,x1t(0)=1,x2t(0)=2,z1(0)=z10,z2(0)=alpha},type=numeric,compile=true,'parameters'=[alpha]):
 > infolevel[all]:=0;
 (10)
 > sol('parameters'=[zs[1]]);
 (11)
 > sol(0.2);
 (12)
 > sol('parameters'=[zs[2]]);
 (13)
 > sol(0.2);
 (14)
 > sol('parameters'=[zs[3]]);
 (15)
 > sol(0.2);
 (16)
 > with(plots):
 > odeplot(sol,[t,z1(t)],0..3);
 Warning, cannot evaluate the solution further right of .21144081, probably a singularity
 > odeplot(sol,[t,z2(t)],0..3);
 Warning, cannot evaluate the solution further right of .21144081, probably a singularity
 >

## The procedure form works for this proble...

method =lsode,

see below.

 > restart:
 > ff:= x1 -> evalf(Int(exp(t), t= 0 .. x1, method = _d01ajc) + 1);
 (1)
 >
 > p:=proc(N,x,Y,F) global ff1; local z1; z1:=ff(x); F[1]:=z1; end proc;
 (2)
 > ics:=Array([1]);
 (3)
 > #infolevel[all]:=10;
 > sol:=dsolve(numeric,number=1,procedure=p,initial=ics,procvars=[y(x)],start=0,method=lsode):
 > sol(0.1);
 (4)
 > sol(1);
 (5)
 >
 >

## the only way I know is to use known = ff...

As I mentioned in the other post you are refering to.

The other approach is to use the procedure form of dsolve where I have shown to integrate a right hand side (rhs) which calls a simple function of x*f(x)^2.

p1 can be your procedure where you do the complicated function and then you can call than in the main procedure p. The difficulty with the procedure form is that compile=true may not work (if p1 is not compiled), events/parametric form is a pain (almost impossible).

The attached code shows that you get the same results compared to direct numerical integration.

 > restart:
 > p1:=proc(x,f) x*f^2; end proc;
 (1)
 > p:=proc(N,x,Y,F) global p1; local z1; z1:=p1(x,Y[1]); F[1]:=z1; end proc;
 (2)
 > ics:=Array([1]);
 (3)
 > #infolevel[all]:=10;
 > sol:=dsolve(numeric,number=1,procedure=p,initial=ics,procvars=[y(x)],start=0):
 > sol(1);
 (4)
 > eq2:=diff(y(x),x)=y(x)^2*x,y(0)=1;
 (5)
 > sol2:=dsolve({eq2,y(0)=1},type=numeric):
 > sol2(1);
 (6)
 >

Bogo,

Please see improved code. I replaced the integral procedure with a numerical dsolve solution in compiled parametric form. You can see that the modified version is much faster.

I know of at least one way to improve the speed even further, I will try that soon.

Hope you find this useful.

 > restart:
 > Digits:=15;
 (1)
 > de := diff(Y(x), x)+10^5*(Y(x)^2-YFD(x)^2)/x^2 = 0;
 (2)
 > YFD:=proc (x)  local xi,g; g:=2; if not type(evalf(x),'numeric') then       'procname'(x)    else Re(evalf((1/2)*g*(int((xi^2-x^2)^(1/2)*xi/(exp(xi)+1), xi = x .. infinity))/Pi^2)); end if; end proc;
 (3)
 > IC:=Y(1)=YFD(1);
 (4)
 > #infolevel[all]:=10;
 > sol:=dsolve({de,IC},type=numeric,known='YFD',method=lsode):
 > sol(1);
 (5)

Solution will be evaluated at a later stage

 > #sol(1.00001);
 > #sol(1.0001);

The itnegral equation is converted to a differential equation and solved in parametric and compiled form with respect to x and xi = 100 is taken to be infinity (not changing much beyond x = 45).

 > eq2:=diff(y(xi),xi)=(xi^2-1.*x^2)^(1/2)*xi/(exp(xi)+1.)/Pi^2;
 (6)
 > sol1:=dsolve({eq2,y(x)=0},type=numeric,'parameters'=[x],compile=true,abserr=1e-15,relerr=1e-14):
 Warning, relerr increased from .1e-13 to 5.00001e-013
 > sol1('parameters'=[1]);
 (7)
 > sol1(1);sol1(2);
 (8)
 > sol1(45.00000000000);
 (9)
 > sol1(100);
 (10)
 > YFD2:=proc (x)   global sol1,y; local xi,z1,g; g:=2; if not type(evalf(x),'numeric') then       'procname'(x)    else    sol1('parameters'=[x]); z1:=sol1(100); rhs(z1[2]);#subs(sol1(100),y(xi)):    end if; end proc;
 (11)
 > YFD(1);YFD2(1);
 (12)

One can note that the dsolve solution for the integral is much faster than the integral calucation and provides the same value.

 > de2:=diff(Y(x), x)+10^5*(Y(x)^2-YFD2(x)^2)/x^2 = 0;;
 (13)
 > sol2:=dsolve({de2,IC},type=numeric,known='YFD2',method=lsode):
 > sol2(1);
 (14)

Next time taken is compared for the original dsolve with integral for the procedure and the new dsolve.

 > t11:=time():sol(1.00001);time()-t11;
 (15)
 > t11:=time():sol2(1.00001);time()-t11;
 (16)
 > t11:=time():sol(1.0001);time()-t11;
 (17)
 > t11:=time():sol2(1.0001);time()-t11;
 (18)
 > with(plots):
 > odeplot(sol2,[x,Y(x)],1..2);
 >

Note that in sol2, compile=true can't be used as it is not able to recognize YDF2. If that can be enabled this code will run even faster.