///////////////////////////////////////////////
//
// modelisation d'une corde avec poids 
// positionné en un point
// risque de rupture en fonction 
// de la position du poids
//
///////////////////////////////////////////////////
function y=fctk(x)
    y=ones(x)+rig*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
//
////////////////////////////////////////////////////////////////
// déformation de la corde 
// pour un poids positionné en point quelconque
/////////////////////////////////////////////////////////////////
//
clear all
lambda=20;
rig=0.8 // rig= 0: corde homogène
a=0.3;   // position du poids sur la corde
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;
y=fctk(x);
z=fctp(x-a*ones(x))
figure(1)
clf()
plot(x,y,'thickness',5)
xtitle('tension (non constante) de la corde')
//
figure(2)
clf()
plot(x,z,'thickness',5)
xtitle('force exercée sur la corde, position a=0.3')
//
figure(3)
clf()
plot(x,Ua,'black','thickness',5)
xtitle('deformation de la corde')
//
////////////////////////////////////////////////////////////////
// Risque de rupture de la corde 
// en fonction de la position du poids
/////////////////////////////////////////////////////////////////
//
k12=fctk((1/2)*h)
kN2=fctk((N+1/2)*h)
//
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
//
////////////////////////////////////////////////
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(vala,valR,'red','thickness',5)
xtitle('risque de rupture en fonction de la position du poids')
//

