有关弧长法的内容已经有很多了，程序也有。即便这样，离实用还有很远的路要走。这里再发一个例子。《混凝土结构有限元分析第二版》第245页的例题，书中是手算演示，我这里用python实现。python代码:import mathdef KT(): return 1def Fint(x): if x > 1 : return 1 - 0.15*(x-1) else: return 2*x - x**2f_ext = 1 # 适应书中计算，假设外荷载为1list_lamd = [0.75] # 荷载步只有一步tol = 2e-2 # 收敛参数，这个值主要是适应书中 5 步迭代conv = 2e11max_iter = 200 # 单个荷载步最大迭代步数niter = 0last_sum_lamd = 0.75 # 上一个收敛点，以下过程跟他没有关系？last_fe = last_sum_lamd * f_ext # 上一个收敛点对应的外荷载 1*0.75=0.75last_sum_x = 0.5 # 上一个收敛点print('niter x lamda conv')for delta_lamd in list_lamd: delta_f = f_ext * delta_lamd # 1*0.75 = 0.75 lamd = 1 Kt = KT() # 常刚度迭代，同一荷载步下，切线刚度矩阵相同 sum_x = 0 while ( conv > tol and niter < max_iter ): niter += 1 if niter == 1: x1 = delta_f / Kt sum_x = last_sum_x + x1 f_int = Fint(sum_x) sum_fe = last_fe + lamd * delta_f # f_resid = sum_fe - f_int else: delta_x1 = delta_f / Kt delta_x2 = f_resid / Kt delta_lamd = (x1 * delta_x2)/(x1 * delta_x1 + lamd * delta_f**2) delta_lamd = -1 * delta_lamd lamd = lamd + delta_lamd sum_fe = last_fe + lamd * delta_f # dx = delta_lamd * delta_x1 + delta_x2 x2 = x1 + dx sum_x = last_sum_x + x2 f_int = Fint(sum_x) sum_fe = last_fe + lamd * delta_f # f_resid = sum_fe - f_int x1 = x2 conv = math.sqrt( f_resid**2 ) print(format(niter, '>3d'), format(sum_x, '>15.5f'), format(lamd, '>16.6f'), format(f_resid, '>18.6e') ) last_sum_lamd = sum_fe last_sum_x = sum_xprint(last_sum_lamd)print(last_sum_x) 迭代过程如图所示代码：import numpy as npimport matplotlib.pyplot as pltx1 = np.linspace(0,1,10)y1 = 2*x1-x1**2x2 = np.linspace(1,3,10)y2 = 1 - 0.15*(x2-1)# Figure 并指定大小plt.figure(num=3,figsize=(8,5))# 绘制F图像，设置 color 为 red，线宽度是 2，plt.plot(x1, y1, color='red',linewidth=2.0)plt.plot(x2, y2, color='red',linewidth=2.0)# 四次迭代的数据iter1 = [0.5, 1.25]iter1_0 = [0.75, 1.5]iter2 = [0.5, 1.519]iter2_0 = [0.75, 1.231]iter3 = [0.5, 1.618]iter3_0 = [0.75, 1.023]iter4 = [0.5, 1.64]iter4_0 = [0.75, 0.929]iter5 = [0.5, 1.643]iter5_0 = [0.75, 0.909]plt.plot(iter1, iter1_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)plt.plot(iter2, iter2_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)plt.plot(iter3, iter3_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)plt.plot(iter4, iter4_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)plt.plot(iter5, iter5_0, label="M1", linewidth = 2, color = "green",marker='o', markersize=7)# 设置 x，y 轴的范围以及 label 标注plt.xlim(0,3)plt.ylim(0,1.7)plt.xlabel('x',fontsize = 18)plt.ylabel('F',fontsize = 18)# 显示图像plt.show() 对于荷载因子的理解有些混乱，还得继续查文献。弧长法的Python实现非线性| 弧长法算例非线性|弧长法改进来源：数值分析与有限元编程
