///////////////////////////////////////////////
//
// modelisation d'une corde avec poids positionné en un point
// optimisation de la position du poids
//
///////////////////////////////////////////////////
function y=fctk(x)
    y=ones(x)+0.8*sin(10*x).*cos(5*x)
endfunction
//////////////////////////////////////////////////
function y=fctp(x)
    y=-lambda*exp(-(lambda*x).*(lambda*x))
endfunction
/////////////////////////////////////////////////
function Fa=fctFa(a)
Fa=zeros(N,1);
for i=1:N
   Fa(i)=fctp(i*h-a); 
end
endfunction
////////////////////////////////////////////////////////////////
// forme de la corde pour un poids positionné en point quelconque
/////////////////////////////////////////////////////////////////
lambda=20;a=0.3;N=80;h=1/(N+1);
//
M=zeros(N,N);
    for i=1:N
       M(i,i)=(fctk((i+1/2)*h)+fctk((i-1/2)*h))/(h^2)
    end
    for i=1:N-1
       M(i,i+1)=-fctk((i+1/2)*h)/(h^2)
    end
    for i=2:N
       M(i,i-1)=-fctk((i-1/2)*h)/(h^2)
    end    
//
Ua=linsolve(M,-fctFa(a))
Ua=[0;Ua;0];
//
x=0:h:1;
figure(1)
clf()
plot(x,fctk(x),'thickness',5)
xtitle('tension (non constante) de la corde')
//
figure(2)
clf()
plot(x,fctp(x-a*ones(x)),'thickness',5)
xtitle('force exercée sur la corde, position a=0.3')
//
figure(3)
clf()
plot(x,Ua,'thickness',5)
xtitle('position de la corde')
//
/////////////////////////////////////////////////
// optimisation de la position du poids
/////////////////////////////////////////////////
function y=pprime(x)
    y=2*lambda^3*x.*exp(-(lambda*x).*(lambda*x))
endfunction
/////////////////////////////////////////////////
k12=fctk((1/2)*h);kN2=fctk((N+1/2)*h);
//
function Ga=fctGa(a)
Ga=zeros(N,1);
for i=1:N
   Ga(i)=-pprime(i*h-a); 
end
endfunction
/////////////////////////////////////////////////
function y=fctR(a)
//
Ua=linsolve(M,-fctFa(a))
y=1/2*k12^2*Ua(1)^2/(h^2)+1/2*kN2^2*Ua(N)^2/(h^2) 
//
endfunction
/////////////////////////////////////////////////
function y=RNprime(a)
//
Va=linsolve(M,-fctGa(a))
Ua=linsolve(M,-fctFa(a))
//
y=k12^2*Ua(1)*Va(1)/(h^2)+kN2^2*Ua(N)*Va(N)/(h^2) 
//
endfunction
//////////////////////////////////////////////
rho=0.03;a0=0.3;RPr=RNprime(a0);
a=a0;niter=0
//
while abs(RPr)>1E-6
  a=a-rho*RPr;niter=niter+1;
  RPr=RNprime(a);  
end
disp('position optimale');disp(a)
disp('nombre d''iterations');disp(niter)
//
////////////////////////////////////////////////
vala=0.1:0.01:0.9;
valR=zeros(vala)
for i=1:length(vala)
   valR(i)=fctR(vala(i))  
end
//
Ua=linsolve(M,-fctFa(a))
Ua=[0;Ua;0];
//
figure(4)
clf()
plot(x,Ua,'thickness',5)
xtitle('forme de la corde (position de poids optimale)')
figure(5)
clf()
plot(vala,valR,'thickness',5)
xtitle('fonction R à minimiser')
//


