MATLAB上机实验――采用龙贝格积分计算轴韧性模量

1、数值分析应用案例: 采用龙贝格积分计算轴韧性模量某杆在轴向负载的作用下会发生变形，其应力应变曲线如图所示。曲线下方从应力为0的点到破裂点的面积称为材料的韧性模数。它提供了一种方法，可以测量出要给定单位体积的材料施加多大的能量才能导致材料破裂。因此，它代表着材料承受冲击负载的能力。对于如图所示的应力应变曲线，利用数值积分计算韧性模数。给定点数据表e0.020.050.100.150.200.25s40.037.543。052.060.055.0（a）轴向负载的杆（b）相应的应力应变曲线，其中应力的单位是千镑每平方英寸（103lb/in2）而应变是无量纲量解：采用数值积分计算韧性模数M=00.25

2、f(e)de首先，需要通过实际测量所得的离散数据点对曲线进行多项式拟合。e0.020.050.100.150.200.25s40.037.543。052.060.055.0采用牛顿差值进行函数拟合：（参考教材应用数值分析P173例5-4）>> x = 0.02,0.05,0.10,0.15,0.20,0.25;y = 40,37.5,43,52,60,55;xx = 0,0.03,0.08,0.13,0.18,0.22;yy,c,p = newtonintep(x,y,xx);>> disp(poly2str(p,'x'); -391428.2175 x

3、5 + 243571.1631 x4 - 64178.4962 x3 + 9133.9186 x2 - 517.8208 x + 47.1786>> 牛顿插值程序(参考教材应用数值分析P172算法5-2 Newton插值)function yy,c,p = newtonintep(x,y,xx)%x,y插值原始数据，要求x按照升序排列，x,y长度一致%xx欲插值点，可以为标量或向量%yy插值结果%c可选，是牛顿差值系数%p可选，插值多项式系数%第1步，变量初始化n = length(x);Q = zeros(n,n);%插商表p = zeros(1,n);Q(:,1) = y(:)

4、;%第一列为函数值%第2步，构造插商表for I = 1:n-1for j = 1:i Q(i+1,j+1) = (Q(i+1,j)-Q(I,j)/(x(i+1)-x(i-j+1);endend%第3步，利用秦九昭算法计算牛顿差值结果yy = Q(n,n);for I = n-1:-1:1yy = yy.*(xx-x(i)+Q(I,i);end%第4步，输出if nargout>1c = Q;if nargout>2 s = zeros(1,n); p = zeros(1,n); for I = 1:n t = 1; for j = 1:i-1 t = conv(t,1,-x(j)

5、; end s(1,n-i+1:n) = t*Q(I,i); p = p+s; endendend得到拟合的多项式函数：f(e)=-391428.21e5 + 243571.16e4 64178.49e3 + 9133.92e2 -517.82e + 47.18函数图像为：>>vel = (x) -391428.2175*x5+243571.1631*x4-64178.4962*x3+9133.9186*x2-517.8208*x+47.1786 %创建匿名函数velvel = (x)-391428.2175*x5+243571.1631*x4-64178.4962*x3+9133

6、.9186*x2-517.8208*x+47.1786>> fplot(vel,0,0.25)%通过内置函数fplot（函数,区间）绘制函数图形>>进行积分计算：M=00.25(-391428.21e5+243571.16e464178.49e3+9133.92e2-517.82e+47.18)de采用龙贝格积分公式积分：>>vel=(x) -391428.2175*x5+243571.1631*x4-64178.4962*x3+9133.9186*x2-517.8208*x+47.1786;>> q,ea,iter=romberg(vel,0,

7、0.25)q = 12.156167862141933ea = 0iter = 3>>龙贝格积分程序（参考工程与科学数值方法的MATLAB实现，美Steven C.Chapra 著.P450算例程序）function q,ea,iter=romberg(func,a,b,es,maxit,varargin)%func被积分函数，a,b积分上下限，es精度，默认为0.000001%q积分估计值，ea误差，iter迭代步数%变量初始化n = 1;I(1,1) = trap(func,a,b,n,varargin:);iter = 0;maxit = 50;es = 0.000001;w

8、hile iter<maxit iter = iter+1; n = 2iter; I(iter+1,1) = trap(func,a,b,n,varargin:); for k = 2:iter+1 j = 2+iter-k; I(j,k) = (4(k-1)*I(j+1,k-1)-I(j,k-1)/(4(k-1)-1); end ea = abs(I(1,iter+1)-I(2,iter)/I(1,iter+1)*100; if ea<=es,break; endendq = I(1,iter+1);复合梯形法则程序（参考工程与科学数值方法的MATLAB实现，美Steven C.Chapra 著.P421算例程序）function I = trap(func,a,b,n,varargin)%func被积分函数，a,b积分上下限，n迭代步数%h积分步长x = a;h = (a-b)/n;s = func