I am a little stuck with my optimization assignment. It's, no surprise, the catenary but the approach we are looking for is a little different than whatever I could track down on different forums across Google. Or maybe I just don't get something and can't incorporate it correctly into my code.
Background
The first part was solving the unconstrained problem (soap film minimizing area), and the second part is the fixed length chain case. I am not able to put the whole script in here as it involves a few different functions, and also I need help with only this one specific idea and not sorting out the whole project, so I will try to explain this as well as possible based on only part of the code - based on threads I've encountered in Matlab forum already (https://se.mathworks.com/matlabcentral/answers/531298-finding-minimum-maximum-of-function-using-lagrange-multipliers and https://se.mathworks.com/matlabcentral/answers/436044-is-there-anyone-to-help-me-for-lagrange-multipliers-method) I have this violent presumption that I am VERY close to figuring that out, I'm just missing something important.
Intro
Let's start from the unconstrained input so it's easier to point out the differences:
function y=F_Discrete(U, N)
%Input
X = 3; %Separation between rings
delX = X/(N-1); %Segments definition
U(1) = 2; %Start height = first ring radius
U(N) = 3; %End height = second ring radius
f_1 = (U(1)*sqrt(1+((U(2)-U(1))/(delX))^2)); %First sum component
f_N = (U(N)*sqrt(1+((U(N)-U(N-1))/(delX))^2)); %Last sum component
f = 0; %Rest of sum preallocation
for i=2:(N-1)
f = f + (U(i)*sqrt(1+((U(i+1)-U(i-1))/(2*delX))^2)); %Second to (N-1)th component
end
%Sum
y = delX * (1/2) * (f_1 + 2*f + f_N);
Nothing fancy here, it's basically taking the Area integral and using the trapezoidal rule on it:
I even have a slide from my report to show steps that were behind this one. It works perfectly fine:
*Problem *
Now is the part I am stuck with. The new code I came up with looks like that:
function y=F_Discrete(U, N)
%Input
X = 4; %Separation
delX = X/(N-1); %Segments definition
U(1) = 2; %Start height
U(N) = 2; %End height
L = 6; %Cable length
%Energy
f_1 = (U(1)*sqrt(1+((U(2)-U(1))/(delX))^2)); %First sum component
f_N = (U(N)*sqrt(1+((U(N)-U(N-1))/(delX))^2)); %Last sum component
f = 0; %Rest of sum preallocation
for i=2:(N-1)
f = f + (U(i)*sqrt(1+((U(i+1)-U(i-1))/(2*delX))^2)); %Second to (N-1)th component
end
E = delX * (1/2) * (f_1 + 2*f + f_N);
%Length
g_1 = (sqrt(1+((U(2)-U(1))/(delX))^2)); %First sum component
g_N = (sqrt(1+((U(N)-U(N-1))/(delX))^2)); %Last sum component
g = 0; %Rest of sum preallocation
for i=2:(N-1)
g = g + (sqrt(1+((U(i+1)-U(i-1))/(2*delX))^2)); %Second to (N-1)th component
end
Len = lambda*(delX * (1/2) * (g_1 + 2*g + g_N)) - L == 0
syms lambda
F = E + lambda*Len;
lambda = double(diff(F,lambda));
%Sum
y = E + lambda*Len
% y = E + Len
Don't mind this y = E + Len; in the end, this is only so the whole program runs through (this is only the input for the bigger optimization 'calculator' with different gradient-like methods, in the input N is the number of segments, and U is the value for that specific point that is calculated outside through BFGS or other Conjugated Gradient function), and so I can check how Len or lambda, or whatever I want to check, looks like.
As You can see, I already tried to incorporate the Multiplier in some way here, I also tried a million other things too, deriving it here and there, multiplying different things in different places, but nothing got me even remotely close to getting the values close to the analytical ones. I am sure I just don't understand the concept correctly and my restriction doesn't apply to the value in the way I want it.
Here's where I need help: since I already have the integrated values from the trapezoidal approach, how do I apply the length constrain there? Where to put it? Right now, since Energy value is not related to the constrain anyhow, it doesn't work (lambda = Len from this diff)
My deadline is until the end of the week but I still have a few things to work on there, this is just a bump in the road that totally locks me out! If someone is following the Matlab forum, I already posted it there so sorry for that, I just need to get through this to continue my work. I'm losing my mind a little, please help! :)
EDIT:
@user619894 pointed out (thank You!) that I should describe better WHAT I'm trying to implement, so I'll to give more details on the general idea
I hope that the area code + explanation on the attached picture makes sense, but just to wrap it up in here:
I had $A=\int_0^3y\sqrt(1 + y'^2)dx$ on which I was testing different gradient-like optimization methods, and had to write an input function for (discretize it).
In this case, it was very easy because I simply used the trapezoidal rule and got my result straight away (as on the picture attached).
Now I have exactly the same function (boundaries -2 to 2):
$E = \int_{-2}^2y\sqrt(1 + y'^2)dx$
but with additional twist that is length constraint:
$R = \int_{-2}^2\sqrt(1 + y'^2)dx = 6$
In general, my question is:
How to discretize E in the same way I did with A in the first part, but with including the length constraint?
In the 2nd snipped of code I attached, there are 3 segments:
Trapezoidal rule for integrating E:
%Energy
f_1 = (U(1)*sqrt(1+((U(2)-U(1))/(delX))^2)); %First sum component
f_N = (U(N)*sqrt(1+((U(N)-U(N-1))/(delX))^2)); %Last sum component
f = 0; %Rest of sum preallocation
for i=2:(N-1)
f = f + (U(i)*sqrt(1+((U(i+1)-U(i-1))/(2*delX))^2)); %Second to (N-1)th component
end
E = delX * (1/2) * (f_1 + 2*f + f_N);
Trapezoidal rule for integrating R (length constraint):
%Length
g_1 = (sqrt(1+((U(2)-U(1))/(delX))^2)); %First sum component
g_N = (sqrt(1+((U(N)-U(N-1))/(delX))^2)); %Last sum component
g = 0; %Rest of sum preallocation
for i=2:(N-1)
g = g + (sqrt(1+((U(i+1)-U(i-1))/(2*delX))^2)); %Second to (N-1)th component
end
Len = lambda*(delX * (1/2) * (g_1 + 2*g + g_N)) - L == 0
where I tried to include the LaGrange multiplier
and calculating the multiplier value, that doesn't work as it should:
syms lambda
F = E + lambda*Len;
lambda = double(diff(F,lambda));
I already see that it's wrong: E has nothing to do with lambda at any point so it's not 'constrained' per se.
I hope this makes it a little clearer but please, let me know if You need any more details! I'll try to explain the thing as long as it's necessary, I am aware I may not be the most precise because I don't understand the concept as well as I would like to
EDIT2: I can see a few people viewing the post, hopefully someone is going to join - until then I will just make use of this place to brainstorm a little and update it with new progress!
First, my most recent change of approach. I decided to take a step back and try to tackle it as it should be from the beginning:
if I take under consideration what I wrote before:
$E = \int_{-2}^2y\sqrt(1 + y'^2)dx$
constrained with:
$R = \int_{-2}^2\sqrt(1 + y'^2)dx = 6$
the Lagrange multiplier form is going to look like this:
$\mathcal{L}(x,y,\lambda) = \int_{-2}^2y\sqrt(1 + y'^2)dx + \lambda(\int_{-2}^2\sqrt(1 + y'^2)dx - 6)$
Here is the part I am completly not sure of:
$\frac{\partial{\mathcal{L}}}{\partial{x}} = 0$
$\frac{\partial{\mathcal{L}}}{\partial{y}} = \int_{-2}^2\sqrt(1 + y'^2)dx$
$\frac{\partial{\mathcal{L}}}{\partial{\lambda}} = \int_{-2}^2\sqrt(1 + y'^2)dx - 6$
Because I have no idea how to calculate the partial derivative after y when there is $(\frac{dy}{dx})^2$ inside the definite integral under the square root... The system obviously looks incorrect as there is nothing to solve for so I am doing something wrong, but I have no idea what.
On the positive note, out of these equations, the integral $\int_{-2}^2\sqrt(1 + y'^2)dx$ I have straight away from the trapezoidal rule (Len in code, just without 'lambda' and '-L' mentioned).
My second idea was to somehow incorporate the lambda and -L into the trapezoidal rule summation (in the code, modify g_1, g and g_N) but it didn't work either. It was promising, but it seems to give completely different results depending on the choice of N.
If anyone is having any ideas - I am super happy to hear anything, I'm on the verge of a mental meltdown at this point :)