Brachistochrone in Freefem

General Discussion
1

I have an interest in simple constrained and unconstrained optimization problems. The path of fastest descent between two points for a bead sliding on a wire in a gravity field can be solved in Freefem. I have attached my edp program which uses a 1D mesh with P1 elements to describe this path. The non-linear part of the problem is a separate factor in the integrand and is stored in an fespace function ‘h’, then iterated on using intermediate values of the path solution from the finite element solver. The optimal path (a cycloid) is converged on in ~10 iterations.

Is there an alternative way to get this same result in Freefem?
mybrachi1D.edp (2.4 KB)

2

I have try other method

argmin J(u) = int1d(Th)( sqrt((1+dx(u)^2)/(u)) ) ; 
with BC u(0)=0, u(L)=h

The a fix point converge very slowly, and the Newton give non invertible matrix.

But I go very fast so I can make mistake my script

The numerical difficulty is the solution have infinite derivative on 0,
so a better think is to do the modélisation with no vertical gravity ie. oX is the slope

macro f(u) ((1+dx(u)^2)/(u)) //
macro df(u,v) (2.*dx(u)*dx(v)/u -  v* (1+dx(u)^2)/u^2)//
macro ddf(u,v,w) (2.*dx(w)*dx(v)/u - (2.*dx(u)*dx(v)*w/u^2 -  2*v*dx(w)*dx(u)^2)/u^2 +  2*v*w*(1+dx(u)^2)/u^3)//
macro dfLb(u,up,v) (2.*dx(u)*dx(v)/up)//  Linearisation of df bilinear part 
macro dfLl(up,v)  (-  v* (1+dx(up)^2)/up^2)//  Linearisation of df linear part 
macro dj(u,v) (0.5*df(u,v)/sqrt(f(u))) //  
macro ddj(u,v,w) (0.5*ddf(u,v,w)/sqrt(f(u)) - 0.5* df(u,v)*df(u,w)/f(u)^1.5)  //
macro DJ(u,v) int1d(Th)(dj(u,v))//
macro D2J(u,v,w) int1d(Th)(ddj(u,v,w))//
// initialisation ..
load "ff-ipopt"
load "msh3"
int L=6; // L is the length of the 1D domain (unstretched bar)
real g=9.8; 
meshL Th=segment(60,[6*x^3,0,0]);

fespace Vh(Th,P1);             // Function space Vh with domain Th

Vh u,v,up,u0;
real h = 1;

u =0; 
solve Pb0(u,v)= int1d(Th)(dx(u)*dx(v))+ on(1,u=0.) +on(2,u=h);
u0[]=u[];
meshL Thmv0=movemesh(Th,[x,-u0]);

The loop of Fixed point algo

for (int i=0; i<4; ++i)
{
	up[]=u[]; // previous value ...
	solve PB1(u,v) = int1d(Th)( dfLb(u,up,v)/ sqrt(f(up)) ) + int1d(Th)(dfLl(up,v)/ sqrt(f(up)))
	               +  on(1,u=0.00) +on(2,u=h);
	up[] -= u[];
	real err = up[].linfty;
	meshL Thmv=movemesh(Th,[x,-u]);
	cout <<i <<  " err="<< err << " cout  " << int1d(Th)( sqrt(f(u)) )<<  endl;
	plot(Thmv,Thmv0,dim=2,cmm=" err= "+err);
	if(err<1e-6) break;
}

the loop for Newton method are

for (int i=0; i<100; ++i)
{
	up[]=u[]; // previous value
	solve Tangent(u,v)= int1d(Th)(ddj(up,u,v)) 
	+ int1d(Th)(-ddj(up,up,v) 
	            + dj(up,v))  
	 +  on(1,u=0.00) +on(2,u=h);
 	up[] -= u[];
 	real err = up[].linfty;
 	meshL Thmv=movemesh(Th,[x,-u]);
 	cout <<i <<  " err="<< err << " cout  " << int1d(Th)( sqrt(f(u)) )<<  endl;
 	plot(Thmv,Thmv0,dim=2,cmm=" err= "+err);
 	if(err<1e-6) break;
	
}