数值计算方法.docVIP

第二章作业 冬WEEK2 陈杨 3110100040 目录 1，问题叙述 1 2，问题分析 2 3，问题求解 2 （1）二分法： 2 （2）试位法 4 （3）不动点迭代法： 5 （4）Newton-Raphson法： 7 （5）割线法: 9 （6）不同方法的比较： 10 1，问题叙述 该梁的挠度曲线方程是: 其中，w0=2.5kN/cm， E=50000kN/cm2为弹性模量， I=30000cm4为转动惯量， L=600cm为杆长。请计算出该杆的最大挠度值。 （提示：杆最大挠度值在满足的x处达到） 2，问题分析 3，问题求解二分法： 观察图可知，该方程有2个解。x=600可能是其中的一个解。带入验证得f(x)=0.说明x=600确实为方程的一个解。根据挠度的物理含义，x=600时，由于该点被固定，所以弯曲变形时横截面形心沿与轴线垂直方向的线位移 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*x^3-l^4*x)； %计算最大挠度 feval(g,c); disp([c;err;yc;g(c)]); toc; 其中子函数bisect的代码如下： function [c,err,yc]=bisect(f,a,b,delta) %Input - f is the function input as a string f % - a and b are the left and right endpoints % - delta is the tolerance %Output - c is the zero % - yc= f(c) % - err is the error estimate for c ya=feval(f,a); yb=feval(f,b); if ya*yb0 error(invalid input); else max1=1+round((log(b-a)-log(delta))/log(2)); for k=1:max1 c=(a+b)/2; yc=feval(f,c); if yc==0 a=c; b=c; elseif yb*yc0 b=c; yb=yc; else a=c; ya=yc; end if (b-a)delta, break,end end end c=(a+b)/2; err=abs(b-a); yc=feval(f,c); end （2）试位法2.683281572999704e+002 误差err 0.000000000000455e+002 f(x)的值0.000034027099609e+002 挠度计算式g 0.005151900620160e+002 计算时间(s) 0.003562 迭代次数k 42 与二分法不同，每次区间的划分点不在中点，而是两区间函数值点连线与x轴的交点，递归地做这样的划分求解最终的结果。Matlab程序如下： clear; clc; x_l = 200; % the initial value of the left bound of x is 0 x_r = 300; % the initial value of the right bound of x is 1 x_c = x_r; % x_c is the cross point of the two bound-line and the x-axis x_t = 0; % x_t is the target ROOT of the eqation e_a = 5e-11; % the tolerance i = 0; % iteration tag l=600; f=@(x)5*x.^4-6*l^2*x.^2+l.^4; %define the function to be solved while (abs(x_t-x_c)e_a) x_t=x_c; %