%% Hinweis: % Man könnte auch bei geschickter Wahl der Startwerte ein % Quasi-Newton-Verfahren verwenden um schnellere Konvergenz zu haben. Man % könnte auch noch die Hesse-Matrix per Hand ausrechnen, um ein % Newton-Verfahren zu verwenden. Dies ist allerdings recht aufwendig. function [Tu,u] = Aufgabe2(P1,P2) % Wähle entweder P1 und P2 (dann müss der plot der exakten Lösung angepasst % werden), oder gebe keine Punkte vor, dann wird automatisch P1 = [1,4] und % P2 = [3.9,0] gewählt dbstop if error % Wähle Punkte, falls nicht spezifiziert if nargin == 0 P1 = [1,4]; P2 = [3.9,0]; end n = 100; % Anzahl äquidistanter Teilintervalle h = (P2(1) - P1(1))/n; % Schrittweite % Startvektor bauen - Verbindungsgerade u0 = zeros(1,n+1); for k = 1:n+1 u0(1,k) = P1(2) - (k-1)*(P1(2)-P2(2))/n; end % Lege Powell-Wolfe Parameter fest gamma = 0.01; eta = 0.9; % Toleranz für Grad-Verfahren tol = 1e-06; % Lösung berechnen [Tu,u] = Grad_PW(u0,gamma,eta,h,tol); % plot plot(P1(1):h:P2(1),u(end,:),'r') hold on % Noch die exakte Lösung: % hierbei wird angenommen, dass P1 = [1,4] und P2 = [3.9,0] ist, sonst % wähle einen anderen Radius und bestimmt damit den neuen P2 (vgl. % Konstruktion der Rollkurve) r = 3; x = r*([0:pi/100:2*pi]-[sin(0:pi/100:2*pi)]); y = r*(-1+[cos(0:pi/100:2*pi)]); plot(x+1,y+4,'b') set(gca,'XLim',[1,4]) legend('Näherung','Exakte Lösung') %keyboard; % um Zugriff auf den workspace zu haben - kann entfernt werden %% Gradient + Powell-Wolfe function [Tu,u] = Grad_PW(u0,gamma,eta,h,tol) % Verwende einfach das Gradientenverfahren mit Powell-Wolfe Schrittweite. % Die Funktionalssauswertung und der Gradient sind als separate Funktionen % gegeben. % ersten Schritt speichern % die Randwerte bleiben fest, vgl. Auswertung vom Gradient -> dieser ist an % den Randknoten stets 0, sodass der dortige Funktionswert nie verändert % wird und damit als Randwert konstant gehalten wird. i = 1; u(i,:) = u0; Tu(i,1) = Brachistochrone_eval_T(u(i,:),h); % starte Iteration while (true) direction = -Brachistochrone_eval_gradT(u(i,:),h); % norm(direction); % optional output if norm(direction) < tol % Abbruchskriterium sehr soft gewählt, da es sonst sehr lnge braucht, vgl Hinweis oben und der Fakt, dass man ein gewisses Diskretisierungslevel für sinnvolle Ergebnisse braucht. break end step_width = Powell_Wolfe(u(i,:),direction,gamma,eta,h); % Den Schritt vollziehen und speichern i = i+1; u(i,:) = u(i-1,:) + step_width*direction'; Tu(i,1) = Brachistochrone_eval_T(u(i,:),h); end %% Powell-Wolfe function s = Powell_Wolfe(u,direction,gamma,eta,h) % Realisierung der Powell-Wolfe-Schrittweitenregel % s_minus halbieren, bis es die Armijo-Bedingung erfüllt s_minus = 2; while (true) s_minus = s_minus/2; if Brachistochrone_eval_T(u+s_minus*direction',h) - Brachistochrone_eval_T(u,h) <= gamma*s_minus*Brachistochrone_eval_gradT(u,h)'*direction; % Armijo-Bedingung break end end % s_plus verdoppeln, bis es die Armijo-Bedingung verletzt s_plus = s_minus; while (true) s_plus = 2*s_plus; if Brachistochrone_eval_T(u+s_plus*direction',h) - Brachistochrone_eval_T(u,h) > gamma*s_plus*Brachistochrone_eval_gradT(u,h)'*direction; % negierte Armijo-Bedingung break end end % Intervallhalbierung bis die Powell-Wolfe-Bedingung erfüllt ist while Brachistochrone_eval_gradT(u+s_minus*direction',h)'*direction < eta*Brachistochrone_eval_gradT(u,h)'*direction s_star = (s_plus + s_minus)/2; % Intervallmitte % Falls Intervallmitte Armijo erfüllt, rechte Intervallhälfte wählen if Brachistochrone_eval_T(u+s_star*direction',h) - Brachistochrone_eval_T(u,h) <= gamma*s_star*Brachistochrone_eval_gradT(u,h)'*direction; % Armijo-Bedingung s_minus = s_star; else % andernfalls linke Intervallhälfte wählen s_plus = s_star; end end s = s_minus; %% evaluate function function [z] = Brachistochrone_eval_T(u,h) % Auswertung des diskreten Funktionals T_h n = length(u); % Anzahl Knoten NOI = n-1; % number of intervalls val = 0; g = 9.81; % Wir nehmen wie immer g = 9.81 an % Aufsummierung for i = 1:NOI val = val + sqrt(h^2+(u(i+1)-u(i))^2)/(sqrt(u(1)-u(i))+sqrt(u(1) - u(i+1))); end %z = h*val; % h war die Schrittweite z = sqrt(2/g)*val; %% evaluate gradient function [gradT] = Brachistochrone_eval_gradT(u,h) % Auswertung des Gradienten des diskreten Funktionals T_h n = length(u); % Anzahl Knoten, aber n-2 Freiheitsgrade -> Gradient am Rand = 0 g = 9.81; % Wir nehmen wie immer g = 9.81 an gradT = zeros(n,1); % Gradientenauswertung - Gradient am Rand = 0! for q = 2:n-1 Summand1 = -(u(q+1) - u(q))/(sqrt(u(1)-u(q))+sqrt(u(1)-u(q+1)))*1/sqrt(h^2+(u(q+1)-u(q))^2) + sqrt(h^2+(u(q+1)-u(q))^2)/(sqrt(u(1)-u(q))+sqrt(u(1)-u(q+1)))^2 * 1/(2*sqrt(u(1) - u(q))); Summand2 = (u(q) - u(q-1))/(sqrt(u(1)-u(q-1))+sqrt(u(1)-u(q)))*1/sqrt(h^2+(u(q)-u(q-1))^2) + sqrt(h^2+(u(q)-u(q-1))^2)/(sqrt(u(1)-u(q-1))+sqrt(u(1)-u(q)))^2 * 1/(2*sqrt(u(1) - u(q))); gradT(q) = sqrt(2/g)*(Summand1 + Summand2); end if gradT == real(gradT) % do nothing else warning('gradient is complex!') % give a warning end