#This file gives the solution of the simply-supported beam problem,
#following the strategy of Section 5.2.2.
#
#
phi:=C1*x^2+C2*x*y+C3*y^2+C4*x^3+C5*x^2*y+C6*x*y^2+C7*y^3+C8*x^4+C9*x^3*y+C10*x^2*y^2+C11*x*y^3+C12*y^4+C13*x^5+C14*x^4*y+C15*x^3*y^2+C16*x^2*y^3+C17*x*y^4+C18*y^5;
#
#
sxx1:= diff(phi,y,y); 
syy1:= diff(phi,x,x); 
sxy1:= -diff(phi,x,y); 
#
#The command `unapply(...,x,y) configures the stresses as functions of x,y,
#so that we can find the value at (e.g.) y=b.
#
sxx2:=unapply(sxx1,x,y);
syy2:=unapply(syy1,x,y);
sxy2:=unapply(sxy1,x,y);
#
#We now find the tractions on y=b as
#
t1:=syy2(x,b);
t2:=sxy2(x,b);
#
#and on y=-b
#
t3:=syy2(x,-b);
t4:=sxy2(x,-b);
#
#On x=a, we have
#
t5:=sxx2(a,y);
t6:=sxy2(a,y);
#
#and on x=-a
#
t7:=sxx2(-a,y);
t8:=sxy2(-a,y);
#The stress function is order 5, so the stresses are order 3 in x and y.
#The tractions on y=+b or -b might therefore be polynomials in x of order 3.
#
#We calculate the coefficients of each power of x in these expressions as
#
s1:=coeff(t1,x,3);
s2:=coeff(t1,x,2);
s3:=coeff(t1,x,1);
s4:=coeff(t1,x,0);
s5:=coeff(t2,x,3);
s6:=coeff(t2,x,2);
s7:=coeff(t2,x,1);
s8:=coeff(t2,x,0);
s9:=coeff(t3,x,3);
s10:=coeff(t3,x,2);
s11:=coeff(t3,x,1);
s12:=coeff(t3,x,0);
s13:=coeff(t4,x,3);
s14:=coeff(t4,x,2);
s15:=coeff(t4,x,1);
s16:=coeff(t4,x,0);
#
#The biharmonic equation is 4th order, so applying it to a 5th order polynomial
#generates a first order polynomial.
#
biharm:=diff(phi,x$4)+diff(phi,y$4)+2*diff(phi,x,x,y,y);
#
#which has to be zero for all x and y.
#
#We enforce this by first separating the coefficients of x,y and constants.
#
b1:=coeff(biharm,x,1);
b2:=coeff(biharm,y,1);
b3:=coeff(coeff(biharm,x,0),y,0);
#
#We also calculate the three force resultants on x=a as
#
Fx1:=int(t5, y=-b..b);
Fy1:=int(t6, y=-b..b);
M1:=int(t5*y, y=-b..b);
#
#and on x=-a
#
Fx2:=int(t7, y=-b..b);
Fy2:=int(t8, y=-b..b);
M2:=int(t7*y, y=-b..b);
#
#We now solve for the constants so as to satisfy (i) the strong boundary
#conditions, (ii) the biharmonic equation and (iii) the weak boundary
#conditions.
#
solution:=solve({s1=0,s2=0,s3=0,s4=-p,s5=0,s6=0,s7=0,s8=0,s9=0,s10=0,s11=0,s12=0,s13=0,s14=0,s15=0,s16=0,b1=0,b2=0,b3=0,Fx1=0,M1=0,Fy1=p*a},{C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,C18});
#
#Notice that there are more equations than there are constants. Some of the 
#equations are not linearly independent. However, Maple can handle this if
#there is a solution.
#
#Substitute the solution into the original stress function.
#
phi1:=subs(solution,phi);
#
#and calculate the final stresses as
#
sxx3:=diff(phi1,y,y); 
syy3:=diff(phi1,x,x); 
sxy3:=-diff(phi1,x,y); 
#