% ビームの静的変形（断面積が一定でない）

L = 10;
A = @(x) (2 - (1/10).*x); % 断面積を表す関数
V = @(xi,xj) (2*(xj-xi) - (1/20)*(xj.^2-xi.^2)); % 小区間で切り取られる領域の体積
n = 6; h = L/n;
E = 200;
fext = 5;

xi = h*[0:n-1]';
xj = h*[1:n]';
Vij = V(xi,xj);

e0 = [Vij;0] + [0;Vij];
em = (-1)*[Vij;0];
ep = (-1)*[0;Vij];
K = (E/(h^2))*spdiags([em e0 ep], -1:1, n+1, n+1); % 三重対角行列
fprintf("剛性行列 K\n");
full(K)

f = zeros(n+1,1);
f(n+1) = fext;
fprintf("ベクトル f\n");
f

a = zeros(n+1,1);
a(1) = 1;
fprintf("ベクトル a\n");
a

mat = [ K,  -a; ...
       -a',  0 ];
vec = [ f; 0 ];
sol = mat \ vec;
un = sol(1:n+1);
lambda = sol(n+2);

fprintf("計算結果\n");
un
lambda

fprintf("伸び\n");
un(2:n+1)-un(1:n)
