www.gusucode.com > symbolic 源码程序 matlab案例代码 > symbolic/FunctionalDerivativeUserGuideExample.m
%% Functional Derivatives Tutorial % This example shows how to use functional derivatives in the Symbolic Math Toolbox(TM) using the example of % the wave equation. The wave equation for a string fixed at its ends is solved using functional derivatives. A functional derivative is the % derivative of a functional with respect to the function that the % functional depends on. The Symbolic Math Toolbox(TM) implements % functional derivatives using the <docid:symbolic_ug.bulw228> function. % % Solving the wave equation is one application of functional derivatives. It % describes the motion of waves, from the motion of a string to the % propagation of an electromagnetic wave, and is an important equation in physics. % You can apply the techniques illustrate in this example to applications % in the calculus of variations from solving the Brachistochrone problem to % finding minimal surfaces of soap bubbles. % % Consider a string of length |L| suspended between the two points |x = 0| % and |x = L|. The string has a characteristic density per unit length and a % characteristic tension. Define the length, density, and tension as % constants for later use. For simplicity, set these constants to |1|. Length = 1; Density = 1; Tension = 1; %% % If the string is in motion, the string's kinetic and potential energies % are a function of its displacement from rest |S(x,t)|, which varies with position % |x| and time |t|. If |d| is the density per unit length, the kinetic energy is % % $$T = \int_0^L \frac{d}{2} \biggl(\frac{d}{dt}S(x,t)\biggr)^2 dx.$$ % % The potential energy is % % $$V = \int_0^L \frac{r}{2} \biggl( \frac{d}{dx}S(x,t)\biggr)^2 dx,$$ % % where _r_ is the tension. % % Enter these equations in MATLAB(TM). Since length must be positive, set % this assumption. This assumption allows |simplify| to simplify the % resulting equations into the expected form. syms S(x,t) d r v L assume(L>0) T(x,t) = int(d/2*diff(S,t)^2,x,0,L); V(x,t) = int(r/2*diff(S,x)^2,x,0,L); %% % The action |A| is |T-V|. The Principle of Least Action states that action % is always minimized. Determine the condition for minimum action, by finding the functional % derivative of |A| with respect to |S| using |functionalDerivative| and % equate it to zero. A = T-V; eqn = functionalDerivative(A,S) == 0 %% % Simplify the equation using |simplify|. Convert the equation into its % expected form by substituting for |r/d| with the square of the wave velocity |v|. eqn = simplify(eqn)/r; eqn = subs(eqn,r/d,v^2) %% % Solve the equation using the method of separation of variables. Set % |S(x,t) = U(x)*V(t)| to separate the dependence on position |x| and % time |t|. Separate both sides of the resulting equation using |children|. syms U(x) V(t) eqn2 = subs(eqn,S(x,t),U(x)*V(t)); eqn2 = eqn2/(U(x)*V(t)) tmp = children(eqn2); %% % Both sides of the equation depend on different variables, yet are % equal. This is only possible if each side is a constant. Equate each side % to an arbitrary constant |C| to get two differential equations. syms C eqn3 = tmp(1) == C eqn4 = tmp(2) == C %% % Solve the differential equations using |dsolve| with the condition that % displacement is |0| at |x = 0| and |t = 0|. Simplify the equations to % their expected form using |simplify| with the |Steps| option set to |50|. V(t) = dsolve(eqn3,V(0)==0,t); U(x) = dsolve(eqn4,U(0)==0,x); V(t) = simplify(V(t),'Steps',50) U(x) = simplify(U(x),'Steps',50) %% % Obtain the constants in the equations. p1 = setdiff(symvar(U(x)),sym([C,x])) p2 = setdiff(symvar(V(t)),sym([C,v,t])) %% % The string is fixed at the positions |x = 0| and |x = L|. The condition % |U(0) = 0| already exists. Apply the boundary condition that |U(L) = 0| % and solve for |C|. eqn_bc = U(L) == 0; [solC,param,cond] = solve(eqn_bc,C,'ReturnConditions',true) assume(cond) %% % The solution |S(x,t)| is the product of |U(x)| and |V(t)|. Find the % solution, and substitute the characteristic values of the string into the % solution to obtain the final form of the solution. S(x,t) = U(x)*V(t); S = subs(S,C,solC); S = subs(S,[L v],[Length sqrt(Tension/Density)]); %% % The parameters |p1| and |p2| determine the amplitude of the vibrations. % Set |p1| and |p2| to |1| for simplicity. S = subs(S,[p1 p2],[1 1]); S = simplify(S,'Steps',50) %% % The string has different modes of vibration for different values of |k|. % Plot the first four modes for an arbitrary value of time |t|. Use the % |param| argument returned by |solve| to address parameter |k|. Splot = subs(S,t,0.3); figure(1) hold on grid on tmp = children(S); ymin = double(tmp(3)); for i = 1:4 yplot = subs(Splot,param,i); fplot(yplot,[0 Length]) end ylim([ymin -ymin]) legend('k = 1','k = 2','k = 3','k = 4','Location','best') xlabel('Position (x)') ylabel('Displacement (S)') title('Modes of a string') %% % The wave equation is linear. This means that any linear combination of % the allowed modes is a valid solution to the wave equation. Hence, the % full solution to the wave equation with the given boundary conditions and % initial values is a sum over allowed modes % % $$F(x,t) = \sum_{k=n}^m A_k \sin(\pi k t)\sin(\pi k x),$$ % % where $A_k$ denotes arbitrary constants. % % Use |symsum| to sum the first five modes of the string. On a new figure, % display the resulting waveform at the same instant of time as the % previous waveforms for comparison. figure(2) fplot(subs(1/5*symsum(S,param,1,5),t,0.3),[0 Length]) ylim([ymin -ymin]) grid on xlabel('Position (x)') ylabel('Displacement (S)') title('Summation of first 5 modes') %% % The figure shows that summing modes allows you to model a % qualitatively different waveform. Here, we specified the initial condition is % $S(x,t=0) = 0$ for all $x$. % % You can calculate the values $A_k$ in the equation % $F(x,t) = \sum_{k=n}^m A_k \sin(\pi k t)\sin(\pi k x)$ by specifying a % condition for initial velocity % % $$u_t(x,t=0) = F_t(x,0).$$ % % The appropriate summation of modes can represent any waveform, which is % the same as using the Fourier series to represent the string's motion.