数值计算方法

第二章作业冬WEEK2陈杨3110100040目录TOC\o"1-5"\h\z\o"CurrentDocument"问题叙述1\o"CurrentDocument"问题分析2\o"CurrentDocument"问题求解2二分法：2试位法4不动点迭代法：5Newton-Raphson法：7割线法:9不同方法的比较：101，问题叙述图(a)是载荷为线性分布的均质梁，如图(b)所示:该梁的挠度曲线方程是：WV=0(-X5+2LX3-Lx)120EIL其中，w0=2.5kN/cm,E=50000kN/cm2为弹性模量，I=30000cm4为转动惯量，L=600cm为杆长。请计算出该杆的最大挠度值。(提示：杆最大挠度值在满足空=0的x处达到)dx2，问题分析将原函数求导后，得f=V'=-—二一(L4-6LX2+5X4).则原问题变为求120EILf(x)=0时的x的值，即求L-6Lx2+5x4=0时的x的值。再比较不同x下的y(x)值，取其中的最大值即为问题的解。3，问题求解二分法：在确定方程感兴趣的根值区间后，找到一个两端点函数值异号的子区间，并考察该子区间的中点的函数值，反复利用根的存在性定理，循环减半根估计区间，递归求得确定精度下的根的估计值。在Matlab中先作图来确定根的大致范围，然后再用二分法求精确解。观察图可知，该方程有2个解。x=600可能是其中的一个解。带入验证得f(x)=0.说明x=600确实为方程的一个解。根据挠度的物理含义，x=600时，由于该点被固定，所以弯曲变形时横截面形心沿与轴线垂直方向的线位移即挠度为0，因此这个解应舍去。则由图可知，另一个解x(200,300)。下面用二分法求出x的精确解。Matlab程序如下：formatlong;tic;clc;l=600;f=@(x)5*x.八4-6*l八2*x"2+l"4;%定义要求解的函数a=200;b=300;delta=5e-11;%定义上下界，容限为eps[c,err,yc]=bisect(f,a,b,delta);w0=2.5;e=50000;i=30000;g=@(x)-w0/(120*e*i*l)*(-x"5+2*l八2*xP-l八4*x)；%计算最大挠度feval(g,c);disp([c;err;yc;g(c)]);toc;其中子函数bisect的代码如下：function[c,err,yc]=bisect(f,a,b,delta)%Input-fisthefunctioninputasastring'f'%-aandbaretheleftandrightendpoints%-deltaisthetolerance%Output-cisthezero%-yc=f(c)%-erristheerrorestimateforcya=feval(f,a);yb=feval(f,b);ifya*yb>0error('invalidinput');elsemax1=1+round((log(b-a)-log(delta))/log(2));fork=1:max1c=(a+b)/2;yc=feval(f,c);ifyc==0a=c;b=c;elseifyb*yc>0b=c;yb=yc;elsea=c;ya=yc;endif(b-a)<delta,break,endendendc=(a+b)/2;err=abs(b-a);yc=feval(f,c);end试位法最后运行结果为:项值解c2.683281572999704e+002误差err0.000000000000455e+002f(x)的值yc0.000034027099609e+002挠度计算式g0.005151900620160e+002计算时间Time(s)0.003562迭代次数k42与二分法不同，每次区间的划分点不在中点，而是两区间函数值点连线与x轴的交点，递归地做这样的划分求解最终的结果。Matlab程序如下：clear;clc;x_l=200;%theinitialvalueoftheleftboundofxis0x_r=300;%theinitialvalueoftherightboundofxis1x_c=x_r;%x_cisthecrosspointofthetwobound-lineandthex-axisx_t=0;%x_tisthetargetROOToftheeqatione_a=5e-11;%thetolerancei=0;%iterationtagl=600;f=@(x)5*x.八4-6*l八2*x"2+l"4;%definethefunctiontobesolvedwhile(abs(x_t-x_c)>e_a)x_t=x_c;%updatethetargetvaluewiththecrosspointiff(x_t)==0break;elsex_c=x_r-f(x_r)*(x_l-x_r)/(f(x_l)-f(x_r));%getthecrosspointiff(x_c)*f(x_r)<0x_l=x_c;elsex_r=x_c;endendi=i+1;enderr=x_t-x_c;w0=2.5;e=50000;i=30000;g=@(x)-w0/(120*e*i*l)*(-x"5+2*l八2*xP-l八4*x);%calculatethebiggestdeflectiondisp([i;x_t;err;f(x_t);g(x_t);])运行结果如下:项值解c0.026832815729997e+004误差err0f(x)的值yc0.000000003051758e+004挠度计算式g0.000051519006202e+004计算时间Time(s)0.002216迭代次数i3000不动点迭代法:根据不动点法收敛法则，对于函数f(x)，当lf’(x)|时，不动点法是收敛的。)，作f’(x)图像得:而对于f(x)=5x4一6Lx2+L,(xH0)，作f’(x)图像得:观察到在xe(200,300)时，f'(x)g(-7.8x108,—7x108)，不符合收敛的条件，故此方程不能采用不动点法。下面用Matlab验证用不动点法是发散的：主函数：formatlong;clc;tic;l=600;f=@(x)5火x.八4-6*l八2*x.八2+l.八4;%定义要求解的函数p0=260;tol=5e-11;max1=10000;%定义起始点，容限为5e-11，即10位有效数字，迭代次数设为10000次[k,p,err,P]=fixpt(f,p0,tol,max1);w0=2.5;e=50000;i=30000;g=@(x)-w0/(120火e*i*l)*(-x.八5+2*l八2火x八3-l八4*x);%k^Eaxiz6A0iEfeval(g,p);disp([p;err;f(p);g(p)]);toc;子函数fixtp：function[k,p,err,P]=fixpt(f,p0,tol,max1)%Input-gistheiterationfunction%-p0istheinitialguessforthefixed-point%-tolisthetolerance%-maxiisthemaximumnumberofiterations%Output-kisthenumberofiterationsthatwerecarriedout%-pistheapproximationtothefixed-point%-erristheerrorintheapproximation%-P'containsthesequence(pn}P(1)=p0;fork=2:maxiP(k)=feval(f,P(k-1));err=abs(P(k)-P(k-1));relerr=err/(abs(P(k))+eps);p=P(k);if(err<tol)||(relerr<tol),break;endendifk==maxidisp('maximumnumberofiterationsexceeded')endP=P';解得结果为：p=NaN,err=NaN所以求解是发散的。(4)Newton-Raphson法:迭代法用迭代方程气=&(气一「求解不动点。在应用到对于求方程f(W=0的根上，最朴素的迭代方程构造是：g(x)=f(x)+x当在整个求根区间上满足f任)+1<11时，迭代算法是有效的。为了加快收敛，同时令\g(x)=。・f(x)+x1g'(x)或f'(x)+1=0这时可以最快速度收敛，可以发现，上面式子稍变形就可以得到：f'(x)即得到下面要使用的牛顿方法:f(x)f'(x)由此可见这一类迭代方法的“收敛性”意义上的最优进化即是牛顿方法。对于其他构造迭代方程的方法这里就不再讨论。迭代方程是：_f3）x—x—^―1—'（n—】）采用绝对预定误差E=5-11，使得结果小数点后有10位有效数字。a,dMatlab程序如下：主函数：formatlong;tic;clc;l=600;f=@(x)5火x.八4-6*l八2*x.八2+l.八4;%定义要求解的函数df=@(x)20火x.八3-4320000火x;%定义求解函数的导数p0=200;epsilon=5e-11;delta=5e-11;max1=10000;%定义起始点，容限为5e-11即10位有效数字，最大迭代次数10000次[p0,err,k,y]=newton(f,df,p0,delta,epsilon,max1);w0=2.5;e=50000;i=30000;g=@(x)-w0/(120火e*i*l)*(-x.八5+2*l八2火x八3-l八4*x);%计算最大挠度feval(g,p0);disp([k;p0;err;f(p0);g(p0)]);toc;子函数newton：function[p0,err,k,y]=newton(f,df,p0,delta,epsilon,max1)%Input-fistheobjectfunctioninputasastring'f'%-dfisthederivativeoffinputasastring'df'%-p0istheinitialapproximationtoazerooff%-deltaisthetoleranceforp0%-epsilonisthetoleranceforthefunctionvaluesy%-max1isthemaximumnumberofiterations%Output-p0istheNewton-Raphsonapproximationtothezero%-erristheerrorestimateforp0%-kisthenumberofiterations%-yisthefunctionvaluef(p0)fork=1:max1p1=p0-feval(f,p0)/feval(df,p0);err=abs(p1-p0);relerr=2火err/(abs(p1)+delta);p0=p1;y=feval(f,p0);if(err<delta)||(relerr<delta)||(abs(y)<epsilon),break,endend运行结果如下:

项值解c2.683281572999748e+002误差err0.000000000000001e+002f(x)的值yc-0.000000457763672e+002挠度计算式g0.005151900620160e+002计算时间Time(s)0.004372迭代次数k4(5)割线法:类似于Newton-Raphson方法：用差分估计斜率:即用x-即用x-xf(x)-f(x)ii-1来代替导数。Matlab程序如下:主函数：formatlong;tic;clc;l=600;f=@(x)5火x.八4-6*l八2*x.八2+l.八4;%定义要求解的函数p0=200;p1=300;epsilon=5e-11;delta=5eT1;max1=10000;%定义起始点，容限为5e-11即10位有效数字，最大迭代次数10000次[p1,err,k,y]=secant(f,p0,p1,delta,epsilon,max1);w0=2.5;e=50000;i=30000;g=@(x)-w0/(120火e*i*l)*(-x.八5+2*l八2火x八3-l八4*x);%计算最大挠度disp([k;p1;err;f(p1);g(p1)]);toc;子函数：function[p1,err,k,y]=secant(f,p0,p1,delta,epsilon,max1)%Input-fistheobjectfunctioninputasastring'f'%-p0andp1aretheinitialapproximationstoazerooff%-deltaisthetoleranceforp1%-epsilonisthetoleranceforthefunctionvaluesy%-max1isthemaximumnumberofiterations%Output-p1isthesecantmethodapproximationtothezero%-erristheerrorestimateforp1%-kisthenumberofiterations%-yisthefunctionvaluef(p1)fork=1:max1p2=p1-feval(f,p1)*(p1-p0)/(feval(f,p1)-feval(f,p0));err=abs(p2-p1);

relerr=2火err/(abs(p2)+delta);p0=p1;p1=p2;y=feval(f,p1);if(err<delta)||(relerr<delta)||(abs(y)<epsilon),break,endend运行结果为:项值解c2.683281572999704e+002误差err0.000000000000455e+002f(x)的值yc0.000034027099609e+002挠度计算式g0.005151900620160e+002计算时