数值大作业.docxVIP

目录大角度单摆的运动周期求解31.单摆运动模型建立31.1 小角度模型32.2 大角度模型32.函数编写52.1mysin函数编写（计算正弦值）52.2 mysqrt函数编写（计算开方值）62.3 f(x)函数编写72.3.1 复合梯形积分法72.3.2 Romberg积分法82.3.3 Gauss积分法（两点）92.4 myspline函数编写（计算三次样条函数）93.运动周期求解104.函数拟合134.1 插值点的选取134.2 插值点的选取134.3 结果分析13总结15大角度单摆的运动周期求解摘要：无阻尼条件下单摆的周期是第一类椭圆积分，在小角度下可以通过取近似来计算周期。但是在大角度下上式不成立，在本文中利用数值积分方法对各个摆角下的周期进行求解，同时亦通过求三次样条插值函数来近似原函数。单摆运动模型建立模型假设：单摆作无阻尼运动摆绳的质量可忽略，且在运动中无形变悬挂的小球可以看做质点参数设计：假设一个无阻尼的理想单摆，初始摆角为，摆长为l，小球的质量为m，重力加速度为g模型求解：根据公式：转动惯量·角加速度 = 力矩，得1.1 小角度模型在小角度（一般认为）的情况下，作近似 则运动方程为这是一个典型的常微分方程，可以解得周期为1.2 大角度模型当初始摆角较大（一般认为）时，上述的近似处理 误差很大，不能使用该方法来求解周期。此时需要精确求解。根据微积分和微分方程的知识，可以解得运动周期的精确表达式为设k为精确周期与近似周期之比，即不妨设则当摆角为时，周期比为为了为之后的计算值提供对比，首先用MATLAB内置函数ellipke计算周期比，将其作为标准值，代码如下所示：theta=1:90; %by degreethm=theta*pi/180; %by radiansk=ellipke(sin(thm/2).^2)*2/pi; %resultxlswrite(elli.xls,[theta;to]) %write results to EXCEL document计算结果如下表所示作出摆角与周期比的函数图像，如下所示。可以看到，当角度很小时，周期比约等于1；当角度增大时，周期比偏离1。可见，当单摆作大角度运动时，近似周期公式是不可靠的。函数编写为了更好地体现数值计算，在本文的函数计算中，均使用四则运算，所以需要编自行写正弦函数mysin和开方函数mysqrt。然后在此基础上，用多种数值积分方法求函数的值。mysin函数编写（计算正弦值）主要思想是利用泰勒公式来求sin(x)的多项式逼近。根据泰勒公式，sin(x)的泰勒展开式如下所示：在本函数中，控制计算值达到小数点后10位，函数代码如下所示：%This function calculate the value of sin(x)%Using Taylors expansionfunction [y] = mysin(x)y = 0;tol = 0.5e-10;i = 0;er = 1;while ertol temp=y; y = y+(-1)^i*x^(2*i+1)/factorial(2*i+1); er = abs(temp-y); i = i+1;end2.2 mysqrt函数编写（计算开方值）%This function calculate the square root of inpout x%Argument x must between zero and one function [y] = mysqrt(x)low = 0;high = 1; %the initial lower and upper boundstol = 0.5e-10; %the principle of 10 significant digitsmid = (high+low)/2; %the estimate of zero point valueer = 1;syms tf(t) = t^2-x;while(ertol) if f(mid)==0 break; elseif f(mid)*f(low)0 high=mid; else low=mid; end temp=mid; %store the last estimated value mid=(high+low)/2; er=abs(temp-mid);endy = mid;2.3 f(x)函数编写令其中x是初始摆角，即最大摆角。则周期比可表示为下面分别用复合梯形积分法、Romberg积分法、Gauss积分法，来构造函数。函数的输入变量是摆角（角度值），输出变量是周期比。2.3.1 复合梯形积分法%You should input the maximum swing angle(by degree) as the