计算固体力学课程作业专 业 固 体 力 学 学 号 1131301009 姓 名 尹亚川作业1：（一）、0=+f H ϕ，其中10=f ，)1(108ϕe H +=（1） 试用直接迭代法，Newton-Raphson 方法，修正Newton-Raphson 方法，拟Newton-Raphson 方法进行求解并进行比较。
（2） 用Euler-Newton 法计算，f 分2级求解：（1）直接迭代法： 0=+f H ϕ（1）)(00ϕH H =（2）于是得近似解)()(101f H -=-ϕ（3）重复这一过程，以第i 次近似解求出第i +1次近似解的迭代公式为)(i i H H ϕ=（4） )()(11f H i i -=-+ϕ（5）直到i i ϕϕϕ-=∆+1（6）变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为5次，-0.9996645=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为4次，-0.9996644=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为2次，-0.9996642=ϕ。
（1） Newton-Raphson 方法 0=+f H ϕ（7） 0)()(≠+=-≡=f H R F i i i i i ϕϕϕψψ（8）ii iT i T K K ⎪⎪⎭⎫ ⎝⎛∂∂≡=ϕψϕ)( （9）)()()(11f H K K i i i T i i T i +=-=∆--ϕψϕ（10） ii iT H K ⎪⎪⎭⎫ ⎝⎛∂∂=⎪⎪⎭⎫ ⎝⎛∂∂=ϕϕϕψ （11）i i i ϕϕϕ∆+=+1（12）当i ϕ∆变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为11次，-0.99966411=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为4次，-0.9996644=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为1次，-0.9996641=ϕ。
（2） 修正的Newton-Raphson 方法将Newton-Raphsom 法迭代公式中的i T K 改用初始矩阵)(00ϕT T K K =，就是修正的Newton-Raphsom 法。
仅第一步迭代需要完全求解一个线性方程组，并将TK 存贮起来，以后的每一步迭代都采用公式)()()(1010f H K K i i T i T i +=-=∆--ϕψϕ（13）当i ϕ∆变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为122244次，-0.986216122244=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为21次，-0.99966321=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为1次，-0.9996641=ϕ。
（3） 拟Newton-Raphson 方法K 的修正要满足一下的拟牛顿方程)()()(111i i i i i K ϕψϕψϕϕ-=-+++（14）对于单变量情况，上式中的1+i K 是导数()i ϕϕϕψ=∂∂的近似表达式，实际上就是割线劲度矩阵。
)()()(100100f H K K i i +=-=∆--ϕψϕ（15） 001ϕϕϕ∆+=（16） 010101011)(ψψϕϕψψϕ--=-∆=-K（17） )()(11111f H K +=∆-ϕϕ （18）ii ii i i i Kψψϕϕψϕ--=∆∆=++-+1111)( （19）当i ϕ∆变得充分小，即近似解收敛时，终止迭代。
取10=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为11次，-0.99966411=ϕ。
取00=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为4次，-0.9996644=ϕ。
取10-=ϕ，令0000005.0<∆ϕ，运用matlab 进行编程求解（代码见附录）。
可得迭代次数为1次，-0.9996641=ϕ。
根据结果可知，在精度取0000005.0<∆ϕ时，Newton 法和拟Newton 法迭代次数基本一致，收敛速度较快，而修正的Newton 法迭代次数较多，收敛速度较慢。
不过，Newton 法和拟Newton 法计算量较大，而修正Newton 法计算量较小。
并且，直接迭代法在解决这种简单问题时迭代次数也较少，收敛速度较快。
若本题不考虑迭代次数，而对精度要求较高，建议采用Newton 法和拟Newton 法；若本题对精度要求不高，主要考虑迭代次数，建议采用Newton 法和拟Newton ；法；若本题对精度和迭代次数要求不高，主要考虑计算量，建议采用修正Newton 法。
Euler-Newton 法在增量步内采用Newton 迭代法。
现以0m δ和m δ分别表示第m 级载荷增量时δ的初值和终值，以m R 表示第m 级增量时的R 的终值，则由式（11）得第m 增量步的迭代公式 10-=m m δδ（20） R R m m λ=（21）()())()(11,11,i m m im T i m m im T i m R K F R K ψλλδ-∆=-=∆---- （22）im i m i m δδδ∆+=+1（23）如果每一增量步内只迭代一次，此时 m m δδ=1（24）m m δδ∆=∆0（25）则对第m 增量步有 ())(011,m m m T m R K ψλδ-∆=∆-- （26）m m m δδδ∆+=-1（27）设00=λ，00=R ，00=δ，5.01=∆λ，5.01=λ，5.02=∆λ，12=λ。
设0000005.0<∆m δ，根据Euler-Newton 法基本原理运用matlab 编程（代码见附录）得001=δ，75.011-=δ，775099.212-=δ，-1.00000022=δ，-0.99966432=δ。
于是可得999664.0-=ϕ。
附录：%直接迭代法clear;y0=1;n=0;for i=1:100;y1=-10/(10*(1+exp(8*y0)));d=y1-y0;y0=y1;if (abs(d)>0.0000005);format long,y0n=n+1;nendend%newton-raphsom法clear;y0=1;n=0;for i=1:100;k=10+10*exp(8*y0)+80*y0*exp(8*y0);f=10*y0+10*y0*exp(8*y0);d=1/k*(-10-f);y0=y0+d;if (abs(d)>0.0000005);format long,y0n=n+1;nendend%修正的newton-raphsom法clear;y0=1;a=1;n=0;for i=1:100000000;k=10+10*exp(8*a)+80*a*exp(8*a);f=10*y0+10*y0*exp(8*y0);d=1/k*(-10-f);y0=y0+d;if (abs(d)>0.0000005);y0n=n+1;nendend%拟newton-raphsom法clear;y0=1;n=0;for i=1:100;f0=10*y0+10*y0*exp(8*y0);k0=10+10*exp(8*y0)+80*y0*exp(8*y0);b=1/k0*(-10-f0);a0=f0+10;y1=y0+b;f1=10*y1+10*y1*exp(8*y1);a1=f1+10;k1=(a1-a0)/b;d=y1-y0;k0=k1;y0=y1;f0=f1;if (abs(d)>0.0000005);y0n=n+1;nendend%Euler_Newton法clear;x0=0;n=0;r=-10;a1=0.5;a2=1;x01=x0;f01=10*(1+exp(8*x01))*x01+10;k01=10+10*exp(8*x0)+80*x0*exp(8*x0);d11=1/k01*(a1*r-f01);x11=x01+d11;x1=x11;x02=x1;f02=10*(1+exp(8*x02))*x02+10;k02=10+10*exp(8*x1)+80*x1*exp(8*x1);d12=1/k02*(a2*r-f01);x12=x02+d12;for i=1:100;k=10+10*exp(8*x12)+80*x12*exp(8*x12);f=10*x12+10*x12*exp(8*x12);d=1/k*(-10-f);x12=x12+d;if (abs(d)>0.0000005);format long,x12dn=n+1;nendend（二）、针对软化问题的求解方法 参考文献：《a local arc-length procedure for strain softening 》 （1）弧长法弧长法的约束方程：21(1,2,3,...)T i l i δδ∆⋅==；其中l ∆为弧长；i δ为现在荷载增量步第i 次迭代的总的增量位移。
i δ的计算式：1()(1,2,3,...)ii j j U i δ∆===∑以外部荷载系数增量∆λ作为未知量，增量位移向量采用Ramm 和Crisfeld 写成：()i F i P U U U ∆∆λ=+⋅1(())(())F i i i U K U R U P λ-=-⋅-⋅ 1(())P i U K U P -=-⋅其中：R 为内部力向量；P 为外部力向量；U 为第i 次迭代总的变形向量；i λ为第i 次迭代总的荷载系数。
i U 和i λ通过第i 次迭代后用下式计算：11(),(1,2,3...)i i ii i iU U U i ∆λλ∆λ-+=+⎧=⎨=+⎩ 由上述方程，增量荷载系数表示成：1T PPlU U ∆∆λ=⋅11()(2,3,4,...)T P i F i TP PU U i U U δ∆λ∆λ-⋅+=-=⋅ 约束方程也改写成：211(1,2,3,...)T l i δδ∆⋅==21(2,3,4,...)T i il i δδ∆-⋅== 1TPPlU U ∆∆λ=⋅2111()(2,3,4,...)T Ti i F i Ti Pl U i U ∆δδ∆λδ----⋅+==⋅（2）局部弧长法May 和Duan 进一步提出局部弧长法，认为在非线性处用相对位移代替1,21,32,1,[...,]n n n ∆δδδδδδδδδ-=----则局部弧长中约束方程为211()()mT e i e e l ∆δ∆δ∆=⋅=∑ 荷载增量表达式为11()()mTP e P ee lUU ∆∆λ∆∆==⋅∑1111()()(2,3,4,...)()()mT P ei F ee i m TP e P ee UU i U U ∆∆δ∆∆λ∆λ∆∆-==⋅+=-=⋅∑∑作业2：开挖荷载的求解方法地基开挖时，需要计算开挖荷载，如下图所示：建立x-y-z 坐标系，z 方向垂直向外，基坑的长为a ，宽为b 。