Merge branch 'master' of ssh://git.code.sf.net/p/maxima/code
[maxima.git] / share / calculus / optvar_1.dem
blob03294763e6e63104198e5fea1f08e0ddab0097d7
1 load('optvar)$
2 depends(y,x)$
3 /*
4 The functional may include derivatives of higher that first
5 order.  For example, including the small-deflection energy
6 of bending, shear, and an elastic foundation, the elastic
7 energy of a nonuniform beam with loading  W(X)  is the
8 integral of the following function: */
9 a(x)*'d(y,x,2)**2 + b(x)*'d(y,x)**2 + c(x)*y**2 + w(x)*y$
10 /* to get the Euler-Lagrange equation: */
11 eq:el(%, y, x);
12 /* If the shear & foundation energy are negligible, as
13 is often the case, we may solve this differential equation,
14 for specific A(X) and  W(X), by two successive applications
15 of ODE2.  For example: */
16 first(eq),d,infeval,a(x)=x, b(x)=0, c(x)=0, w(x)=1$
17 bb:%,'d(y,x,3)='d(dy2,x),'d(y,x,4)='d(dy2,x,2);
19 dependencies(dy2(x))$
20 ev(bb,d,eval);
21 ode2(%, dy2, x);
22 ode2(subst([%k1=%k3,%k2=%k4,dy2='d(y,x,2)],%), y, x);
23 /* Now let's impose the boundary conditions
24 Y='D(Y,X)=0 at X=1, Y=0, 'D(Y,X)=1 at X=2: */
25 ic2(%, x=1, y=0, 'd(y,x)=0);
26 ic2(subst([%k3=%k1,%k4=%k2],%), x=2, y=0, 'd(y,x)=1);
27 /* as another specific beam example: */
28 eq,eval,a(x)=1,b(x)=2,c(x)=1,w(x)=1$
29 first(%),d;
30 /* We cannot treat this as 2 successive second-order
31 equations, but the function DESOLVE in the SHARE file
32 DESOLN, already loaded, is applicable to some linear
33 arbitrary-order constant coefficient equations.  First we
34 must convert the equation so that the dependency of  Y  is
35 explicitly indicated: */
36 convert(%, y, x);
37 /* DESOLVE requires all boundary conditions to be at
38 X=0, and it works best if they are specified before calling
39 the function. However, we may overcome this restriction, as
40 illustrated by the following example, where the beam has
41 zero slope and deflection at both ends: */
42 (atvalue(y(x),x=0,0), atvalue('d(y(x),x),x=0,0),
43 atvalue('d(y(x),x,2),x=0,%k1), atvalue('d(y(x),x,3),x=0,%k2))$
44 desolve(%th(2), y(x));
45 ic2(subst(y(x)=y,%), x=1, y=0, 'd(y,x)=0);