function[D,C]=RitzTorsionD(a,b,n)%Thisfunctioncalculatesthetorsionalrigidityofrectangularsection%basedonRitzmethod,inwhichthetrialfunctionis%(x^2-a^2)*(y^2-b^2)*(C1+C2*x^2*y^2+C3*x^4*y^4+C4*x^6*y^6+...)%Usage:%[C,D]=RitzTorsion(a,b,n)%Input:%a------widthofthesection%b------heightofthesection%n------thenumberofthetrialfunctions%Output:%D------thetorsionalrigidity%C------thecoefficientsofthetrialfunctions%%-------------------------------------------------------------%---------------------------------------------------------------ifnargin3a=1;b=1;n=4;endsymsxyC=zeros(n,1);D=0;F0=(x^2-a^2)*(y^2-b^2);F=sym(zeros(n,1));fori=1:nj=2*(i-1);F(i)=F0*x^j*y^j;endFx=diff(F,x);Fy=diff(F,y);A=zeros(n,n);B=zeros(n,1);fori=1:nforj=1:nA(i,j)=eval(2*int(int(Fx(i)*Fx(j)+Fy(i)*Fy(j),x,-a,a),y,-b,b));endB(i)=eval(4*int(int(F(i),x,-a,a),y,-b,b));endC=A\B;D=sum(C.*B/2);return