Carga crítica en una columna (G 7A)

%Apartado 3
L=9; p=1.2;
R=1:0.1:5;
Pcr=(pi^3*p*R.^4)/(2*L);
plot(R,Pcr)

% Calculamos numéricamente los valores de P_cr. El primer autovalor es el que se corresponde con dicha carga (esto está almacenado en la primera columna de la matriz A)
clear, clf;
x0=0;
L=9;
h=0.1;
N=(L-x0)/h;
x=x0:h:L;
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
R=linspace(1,5,N+1);
for i=1:length(R)
    A(i,:)=eig(0.5*pi*1.2*R(i)^4*L.*K);
end
hold on
Pcr=A(:,1);
L2=9; p2=1.2;
R2=1:0.1:5;
P2=(pi^3*p2*R2.^4)/(2*L2);
figure(1)
plot(R2,P2,'g')
figure(2)
plot(R,Pcr,'-')
xlabel('Radio (m)');
ylabel('Carga critica (kN)')
hold off

%Apartado 8
a=0;b=9;
h=0.001;
N=(b-a)/h;
x=a:h:b
E=1;
I=1/4*1.2*pi*0.8^4;
P=0.8846678;
%introducimos la función -x/EI
f=-x./(E*I);
%indica numero de sumando que tomamos de la serie de fourier
G=[3,5,10,40];
Q=0;
hold on
for i=1:4
    for k=1:G(i)
   
        %introucimos los autovalores
        muk=k*pi;
        %autfunciones
        pk=sin(k*pi*x);
        %calculmo los coef de fourier
        ak=trapz(x,f.*pk)/trapz(x,pk.^2);
        %calculo de la aproximación de fourier
        Q=Q+(ak/(muk+P/E*I))*pk;
    end
    plot(Q,x)
end
hold off