forked from chebfun/examples
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Catenary.m
141 lines (121 loc) · 4.13 KB
/
Catenary.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
%% The catenary problem from the calculus of variations
% Hrothgar and Toby Driscoll, April 2014
%%
% (Chebfun example opt/Catenary.m)
% [Tags: #optimization, #variational calculus, #chebfun2]
FS = 'fontsize'; LW = 'linewidth';
%%
% Suppose we wish to solve a classical problem from the calculus of variations
% of the form
%
% $$ \min_{y(t)} J[y] \equiv \int_{-1}^1 f(t,y,y') dt $$
%
% subject to the boundary conditions
%
% $$ y(-1) = \alpha, \quad y(1) = \beta. $$
%
% Problems like this arise in physics, for example, when $J$ is an energy functional.
%%
% A standard example of a variational problem is the catenary problem, which
% is to determine the shape of a hanging rope. In this case $J$ is an energy
% functional taking the form [1]
%
% $$ J[y] = \int f(t,y,y') dt = \int y \sqrt{1 + (y')^2} dt. $$
%%
% To solve this problem with Chebfun, we begin with the functional and
% boundary conditions. We will choose boundary conditions that make it easy to
% verify the known solution.
f = @(y,yp) y .* sqrt(1 + yp.^2);
J = @(y) sum( f(y, diff(y)) );
alpha = cosh(-1);
beta = cosh(1);
%%
% Since the problem is nonlinear, we will need to solve it iteratively.
% As an initial guess to the iteration we'll choose a simple line, and
% note its functional value.
dom = [-1,1];
y0 = chebfun([alpha, beta]', dom);
startJ = J(y0)
%%
% The simplest way to compute extrema of the functional $J[y]$ is to
% find solutions of the Euler-Lagrange equation
%
% $$ \frac{\partial f}{\partial u} -
% \frac{d}{dx}\frac{\partial f}{\partial (u')} = 0. $$
%
% To do so we must compute the first variations of the function $f$,
%
% $$ \frac{\partial f}{\partial u}, \quad \frac{\partial f}{\partial u'}. $$
%%
% Using only information about first variations, we could implement a steepest
% descent method. However, if can also get our hands on the second variations
% of $f$, then we could implement a more efficient Newton method. The second-
% variation equivalent of the Euler-Lagrange equation, sometimes called
% Jacobi's accessory equation, is
%
% $$ \frac{d}{dx}\Big(\frac{\partial^2 f}{\partial(u')^2}\Big) +
% \Big(\frac{d}{dx}\frac{\partial^2 f}{\partial u \partial u'}-
% \frac{\partial^2 f}{\partial u^2}\Big)u = 0. $$
%%
% To find the first and second variations, we will use Chebfun2.
% First we make a chebfun2 of $f$.
dom2 = [1 2 -2 4];
F = chebfun2(f, dom2);
%%
% Now we can compute the first and second variations of $f(y,y')$ by taking
% derivatives of `F`: each `diff(F,[a b])` is a chebfun2 corresponding to
%
% $$ \frac{\partial^{(a+b)} f}{\partial^a (u) \partial^b (u')}. $$
F1 = diff(F, [1 0]);
F2 = diff(F, [0 1]);
F11 = diff(F, [2 0]);
F12 = diff(F, [1 1]);
F22 = diff(F, [0 2]);
%%
% Now we are able to compute, say $\partial f(y,y')/\partial u$ for a
% particular function $y$ by writing `F1(y,diff(y))` for a chebfun `y`.
%%
% The iteration will take Newton steps by solving the accessory equation. Five
% iterations are enough to get an accurate solution.
y = y0;
pref = cheboppref;
pref.plotting = 'off';
pref.display = 'none';
for k = 1:5
% The first and second variations of f.
yp = diff(y);
f1 = F1(y, yp);
f2 = F2(y, yp);
f11 = F11(y, yp);
f12 = F12(y, yp);
f22 = F22(y, yp);
% Compute the next Newton step by solving the accessory equation.
N = chebop( @(x,u) diff(f22.*diff(u)) + (diff(f12)-f11).*u, dom );
N.lbc = 0;
N.rbc = 0;
u = solvebvp(N, (f1-diff(f2)), pref);
% Increment by the step.
y = y + u;
% Print out the current value of the energy functional.
nextJ = J(y)
end
%%
% The exact solution is a simple $\cosh$. Here is the error in our solution.
y_exact = chebfun('cosh(x)', dom);
norm(y - y_exact)
%%
% And the values of the energy functional:
fprintf('\n final J[y]: %.16f\n', J(y))
fprintf('optimal J[y]: %.16f\n', J(y_exact))
%%
% Finally, a plot of the catenary:
plot(y_exact, 'r--', LW,2), hold on
plot(y, 'k-', LW,1)
axis equal
title('Solution to the catenary problem', FS, 12)
xlabel('x', FS, 14), ylabel('y(x)', FS, 14)
legend('exact','computed','location','north')
%% References
%
% 1. Charles Fox, _An Introduction to the Calculus of Variations_,
% Courier Dover Publications, 1987.