数值计算 实验.doc

数值计算 实验学 院： 建筑工程学院 专业班级： 土木工程硕士1301班 学 号： 2111306021 学生姓名： 李 超 授课老师： 陆建芳 提交日期： 2014.5.31 我是一名土木工程专业的学生，下面的四道题目，背景都是源自于本专业的一些问题。有的是我自己编的，有的是我从书上看到的。1、如图所示刚架，受竖向均布荷载q 的作用。试用位移法求解各节点的转角位移 、 、 、 ， 刚架尺寸 和 线刚度 由图给出 。 建立该刚架的位移法典型方程如下： 12i + 2i - = 0 2i + 20i +4 i + = 0 4i + 20i +2 i - = 0 2 i + 12i + = 0手算的精确解为： = = = = 转角的单位为 （ ）机算算法选择“列主元高斯消去法” , Matlab程序如下:A=input(系数矩阵A：n);b=input(向量B：n);B=A b; n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica0, disp(RA=RB，方程组无解.n) returnendif RA=RB if RA=n fprintf(RA=RB=%d，方程组有唯一解.n,n) X=zeros(n,1); for p=1:n-1 t=find(abs(B(p:end,p)=max(abs(B(p:end,p)+p-1; if abs(B(t,p)=abs(B(p,p) l=B(t,:); B(t,:)=B(p,:); B(p,:)=l; end %列主元判断 for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end end %把方程组系数矩阵A化为同解的上三角矩阵 b=B(1:n,n+1); A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); end %从xn至x1逐个求解上三角方程 else disp(请注意：因为RA=RBn，所以此方程组有无穷多解.) return endenddisp(方程组解为：);X系数矩阵A：12 2 0 0;2 20 4 0;0 4 20 2;0 0 2 12向量B：1/3 -1/4 1/4 -1/3 RA=RB=4，方程组有唯一解.方程组解为：ans = 0.0310 -0.0195 0.0195 -0.0310 精确值与计算值非常接近，误差不到1/100 。 该算法能帮助我们结构工程师求解那些阶数不大但手算复杂的线性方程组，而且精度能满足行业规范的要求（5/100），具有很高的实用性。2、 计算如图所示两层框架结构的自振频率，各层质量分别为：= 60t ，= 50t 。第一层层间侧移刚度，第二层层间侧移刚度 。框架各层的层间刚度系数分别为：则得频率方程如下：将上式展开，得： 手算的精确解为： 17.54 rad / s 40.32 rad / s机算算法选择“割线法” , Matlab程序如下:function x,k,y=gexian(f,x0,x1,tol,n) fprintf (x(%2d)=%10.8fn,0,x0) fprintf (x(%2d)=%10.8fn,1,x1) for k=1:n x=x0 - f(x0)*(x1-x0)/(f(x1)-f(x0); err=abs(x-x1); x0=x1; x1=x; y=f(x); fprintf (x(%2d)=%10.8fn,k+1,x) if (err f=inline(0.00003*x.4-0.058*x.2+15); x,k,=gexian(f,8,10,1e-4,1000)x( 0)=8.00000000x( 1)=10.00000000x( 2)=19.94306288x( 3)=17.36665594x( 4)=17.53755446x( 5)=17.53689426x( 6)=17.53689451x = 17.5369k = 5 x,k,=gexian(f,35,36,1e-4,1000)x( 0)=35.00000000x( 1)=36.00000000x( 2)=43.81280308x( 3)=39.18535863x( 4)=40.06328253x( 5)=40.34434956x( 6)=40.32064789x( 7)=40.32109377x( 8)=40.32109453x = 40.3211k = 7 精确值与计算值非常接近，误差不到1/100 。该算法能帮助我们结构工程师求解某些方程的根，具有很高的实用性。3、 曲杆的位移计算： 如图所示是一等截面的圆弧形曲杆AB,截面为矩形，圆弧AB的圆心角为 ，半径为R . 设均布的竖向荷载 q 沿水平线作用。试求B点的竖向位移 。解：求时，我们在B点加单位竖向荷载，如上图。然后分别求实际荷载和单位荷载作用下的内力。取B点作为坐标原点，任一点C的坐标为 x ，y ，圆心角为 。由虚功原理可得：位移 取作变数，即 ， ，代入上式，得 令 ，精确解经手算，保留四位小数后为： 机算算法选择“复化辛普森”法 , Matlab 程序如下function S_n = comp_simpson(f,a,b,n)h = (b-a)/n;s = 0;for i =0:n-1 s = s + f(a+i*h) + 4*f(a+(i+1/2)*h) + f(a+(i+1)*h);endS_n = h*s/6; f=inline(1/2*(sin(x).3); comp_simpson(f,0,1.3090,3)ans =0.2067即得 精确值与计算值非常接近，误差不到1/100 ，具有很高的实用性。能帮助我们求解一些涉及到积分的工程运算。4、混凝土的应力-应变关系是研究钢筋混凝土构件强度、裂缝和变形理论的重要依据。在一次单调加载试验中获得多组应力应变组合值如下表所示。试用二次曲线来拟合混凝土的 - 曲线。应力510152030应变（）0.20.50.91.22.0算法选择“最小二乘多项式”法 , Matlab 程序如下function p = polyfit(x,y,n)G = zeros(n+1,n+1);for k =0:n for j = 0:n G(k+1,j+1) = sum(x.(k+j); end b(k+1) = sum(x.k.*y);enda = Gb;p = fliplr(a);输入下面的命令： x=0.2 0.5 0.9 1.2 2.0; y=5 10 15 20 30; p=polyfit(x,y,2)p = -0.9326 15.9012 1.9547所以拟合曲线为： = - 0.9326+15.9012+1.9547 （）误差分析：为了更直观地观察曲线拟合的效果，用下面的代码将曲线和原来的数据画在图中，如下截图所示： xi=0:0.1:2.5; yi=polyval(p,xi); plot(x,y,ko, xi,yi,k-)由图可以看出用圆圈表示的数据点紧密地