4609.sagews

# -*- coding: utf-8 -*-
from sympy import *
from sympy.matrices.matrices import *

x = Symbol('x')
P = Symbol('P',positive=True)
L = Symbol('L', positive=True)
EI = Symbol('EI', positive=True)

a = Symbol('a')
b = Symbol('b')

def ele_solve(y,unknown):
 print "Start calculating Pcr for:",y,"with method1, which means using y'' rather than M"
 dy = diff(y,x)
 dy2 = diff(y,x,2) 
 D = integrate((diff(y,x))**2/2 , (x,0,L))
 print "\tDelta = \n\t\t",D.evalf()
 PI = integrate(EI*dy2**2/2,(x,0,L))-P*D
 print "\tPI = \n\t\t",PI.evalf()
 EQs = Matrix([PI]).jacobian(unknown)
 print "\t平衡方程EQs: "
 for e in EQs:
  print "\t\t",e.evalf(),"= 0 "
 EQ = simplify(Matrix(EQs).jacobian(unknown).det())
 print "\t二阶变分为零:"
 print "\t\t",EQ.evalf(),"= 0"
 ans = solve( EQ, P)
 print "\t解得:"
 for p in ans:
  print "\t\tPcr = ",p.evalf()
 model = []
 Pcr = ans
 for p in Pcr:
  bcr = solve( EQs.subs(P,p).subs(a,1)[1:] , unknown[1:])
  model.append([1,bcr])
 return (Pcr,model)

def ele_solve2(y,unknown):
 print "Start calculating Pcr for:",y,"with method2, which means using M rather than y''"
 dy = diff(y,x)
 D = simplify(integrate((dy)**2/2 , (x,0,L)))
 print "\tDelta = \n\t\t",D.evalf()
 
 M = simplify(P*(y.subs(x,L) - y))
 print "\tM = \n\t\t",M.evalf()
 PI = integrate(M**2/EI/2,(x,0,L))-P*D
 print "\tPI = \n\t\t",PI.evalf()
 
 EQs = Matrix([PI]).jacobian(unknown)
 print "\t平衡方程EQs: "
 for e in EQs:
  print "\t\t",e.evalf(),"= 0 "
 
 EQ = Matrix(EQs).jacobian(unknown).det()
 print "\t二阶变分为零:"
 print "\t\t",EQ.evalf(),"= 0"
  
 if len(unknown)>1:
  ans = solve( EQ , P)
  ans.remove(0)
  print "\t解得:"
  for p in ans:
   print "\t\tPcr = ",p.evalf()
 
  model = []
  Pcr = ans
  for p in Pcr:
   bcr = solve( EQs.subs(P,p).subs(unknown[0],1)[1:] , unknown[1:])
   model.append([1,bcr])
  return (Pcr,model)
 else :
  ans = solve(EQs[0], P)
  Pcr = ans
  Pcr.remove(0)
  print "\t解得:"
  for p in Pcr:
   print "\t\tPcr = ",p.evalf()
  return (Pcr,[1])

def result(y,uk,fun):
 (Pcr,Model) = fun(y ,uk)
 Pcr = [ simplify(i) for i in Pcr]
 if 0 in Pcr:
  Pcr.remove(0)
 return min(Pcr)
 
P0 = result(a*(1-cos(pi*x/2.0/L)),[a],ele_solve).evalf()
print "！此解为真实解"
print "\n\n"

def calandanalyse(y):
 P1 = result(y[0],y[1],ele_solve)
 print "\n\n"
 P2 = result(y[0],y[1],ele_solve2)
 print "误差分析：\n"
 print "第一种方法误差:",(abs(P1/P0)-1).evalf()*100,"\t%"
 print "第二种方法误差:",(abs(P2/P0)-1).evalf()*100,"\t%\n"
 print " "*10,"-"*10,"华丽的分割线","-"*10,'\n\n'

calandanalyse([a*x**2,[a]])
calandanalyse([a*(-x**3/L + x**2),[a]])
calandanalyse([a*x**2+b*x**3,[a,b]])
calandanalyse([a*x**2,[a]])

WARNING: Output truncated!

full_output.txt

Start calculating Pcr for: a*(-cos(0.5*pi*x/L) + 1) with method1, which means using y'' rather than M Delta = 0.616850275068085*a**2/L PI = 1.52201704740629*EI*a**2/L**3 - 0.616850275068085*P*a**2/L 平衡方程EQs: 3.04403409481258*EI*a/L**3 - 1.23370055013617*P*a/L = 0 二阶变分为零: 3.04403409481258*EI/L**3 - 1.23370055013617*P/L = 0 解得: Pcr = 2.46740110027234*EI/L**2 ！此解为真实解 Start calculating Pcr for: a*x**2 with method1, which means using y'' rather than M Delta = 0.666666666666667*L**3*a**2 PI = 2.0*EI*L*a**2 - 0.666666666666667*L**3*P*a**2 平衡方程EQs: 4.0*EI*L*a - 1.33333333333333*L**3*P*a = 0 二阶变分为零: 4.0*EI*L - 1.33333333333333*L**3*P = 0 解得: Pcr = 3.0*EI/L**2 Start calculating Pcr for: a*x**2 with method2, which means using M rather than y'' Delta = 0.666666666666667*L**3*a**2 M = L**2*P*a - P*a*x**2 PI = -0.666666666666667*L**3*P*a**2 + 0.266666666666667*L**5*P**2*a**2/EI 平衡方程EQs: -1.33333333333333*L**3*P*a + 0.533333333333333*L**5*P**2*a/EI = 0 二阶变分为零: -1.33333333333333*L**3*P + 0.533333333333333*L**5*P**2/EI = 0 解得: Pcr = 2.5*EI/L**2 误差分析： 第一种方法误差: 21.5854203708053 % 第二种方法误差: 1.32118364233778 % ---------- 华丽的分割线 ---------- Start calculating Pcr for: a*(x**2 - x**3/L) with method1, which means using y'' rather than M Delta = 0.0666666666666667*L**3*a**2 PI = 2.0*EI*L*a**2 - 0.0666666666666667*L**3*P*a**2 平衡方程EQs: 4.0*EI*L*a - 0.133333333333333*L**3*P*a = 0 二阶变分为零: 4.0*EI*L - 0.133333333333333*L**3*P = 0 ... Start calculating Pcr for: a*x**2 + b*x**3 with method2, which means using M rather than y'' Delta = 0.9*L**5*b**2 + 1.5*L**4*a*b + 0.666666666666667*L**3*a**2 M = L**3*P*b + L**2*P*a - P*a*x**2 - P*b*x**3 PI = -0.9*L**5*P*b**2 - 1.5*L**4*P*a*b - 0.666666666666667*L**3*P*a**2 + 0.321428571428571*L**7*P**2*b**2/EI + 0.583333333333333*L**6*P**2*a*b/EI + 0.266666666666667*L**5*P**2*a**2/EI 平衡方程EQs: -1.5*L**4*P*b - 1.33333333333333*L**3*P*a + 0.583333333333333*L**6*P**2*b/EI + 0.533333333333333*L**5*P**2*a/EI = 0 -1.8*L**5*P*b - 1.5*L**4*P*a + 0.642857142857143*L**7*P**2*b/EI + 0.583333333333333*L**6*P**2*a/EI = 0 二阶变分为零: 0.15*L**8*P**2 - 0.0671428571428571*L**10*P**3/EI + 0.00257936507936508*L**12*P**4/EI**2 = 0 解得: Pcr = 2.46804416707517*EI/L**2 Pcr = 23.5627250636941*EI/L**2 误差分析： 第一种方法误差: 0.752232737739857 % 第二种方法误差: 0.0260625158497754 % ---------- 华丽的分割线 ---------- Start calculating Pcr for: a*x**2 with method1, which means using y'' rather than M Delta = 0.666666666666667*L**3*a**2 PI = 2.0*EI*L*a**2 - 0.666666666666667*L**3*P*a**2 平衡方程EQs: 4.0*EI*L*a - 1.33333333333333*L**3*P*a = 0 二阶变分为零: 4.0*EI*L - 1.33333333333333*L**3*P = 0 解得: Pcr = 3.0*EI/L**2 Start calculating Pcr for: a*x**2 with method2, which means using M rather than y'' Delta = 0.666666666666667*L**3*a**2 M = L**2*P*a - P*a*x**2 PI = -0.666666666666667*L**3*P*a**2 + 0.266666666666667*L**5*P**2*a**2/EI 平衡方程EQs: -1.33333333333333*L**3*P*a + 0.533333333333333*L**5*P**2*a/EI = 0 二阶变分为零: -1.33333333333333*L**3*P + 0.533333333333333*L**5*P**2/EI = 0 解得: Pcr = 2.5*EI/L**2 误差分析： 第一种方法误差: 21.5854203708053 % 第二种方法误差: 1.32118364233778 % ---------- 华丽的分割线 ----------