matlab程序设计实践

wR1 =
3.0131e-010
以 作为停止计算的条件，递归5层。
利用龙格贝尔求积公式的roberg.m程序计算I的近似值R=0.3554，
所得的相邻两次迭代的绝对误差为wugu=1.4666e-010
定积分的精确值的近似值Fs=0.3554
fi与R的绝对误差wR=4.4797e-010，已达到所要求的精度1.e-8。
利用Romberg求积公式计算 ，并与分析其误差
1)算法说明
通过Romberg求积公式是一种加速计算积分的方法,在变步长的求积过程中,可以将精度低的梯形值Tn,逐步加工成精度较高的辛卜生值Sn,牛顿-柯茨值Cn与龙贝格值Dn.
设函数 闭区间 上有定义且在n+1个节点 处的函数值 义，则在闭区间 上用二分法每次将n增加一倍，可的Romberg求积公式
apply2.m
x(1)=0;y(1)=0;
x=zeros(1,30000);y=zeros(30000,1);
fori=1:30000
x(1,i+1)=1+y(i,1)-1.4*x(1,i).^2;
y(i+1,1)=0.3*x(1,i);
end
plot(x,y,'b.')
运行结果
3） 流程图
是
否
3、利用数值积分法求
中南大学
MATLAB程序设计实践
菲 儿
材料科学与工程学院
2012年4月11号
1、编程实现以下科学计算算法，并举一例应用之（参考书籍《精通matlab科学计算》，王正林等著，电子工业出版社，2009年）——“牛顿科茨公式”
1)算法说明
在构造矩形公式、梯形公式和辛普森公式时，我们采用简单多项式 代替被积函数 来构造求积公式，因为多项式不但计算方便，而且容易积分。同样的，牛顿-科茨公式即是在 上采用拉格朗日插值多项式 代替被积函数 ，然后求其积分。
0.3476 0.3567000
0.3535 0.3555 0.355400
0.3550 0.3554 0.3554 0.35540
0.3553 0.3554 0.3554 0.3554 0.3554
R =
0.3554
wugu =
1.4666e-010
h =
0.0625
Fs =
0.3554
wR =
4.4797e-010
f=subs(f,'t',x0);
else
f=vpa(f,6);%对逼近值f取六位有效数字
end
end
end
第二个源程序：
>> f=Legendre('1/(2-x)',6)
f=Legendre('1/(2-x)',6，0.5)
③运行结果：
f =
0.0955158*t - 0.0000116697*t*(2.2*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375) - 5.86667*t + 2.93333*t*(7.5*t^2 - 2.5)) - 0.0000565953*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375) + 0.00154825*t*(7.5*t^2 - 2.5) - 0.000275682*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) +0.0305351*t^2 + 0.539128
fi=int(exp((-x.^3)./2)./(sqrt(2*pi)),x,0,1);%用符号计算函数的精确积分值
Fs=double(fi)%转为双精度型
wR=double(abs(fi-R))%计算fi与R的绝对误差
wR1=wR-wugu%计算wR与wugu的绝对误差
运行结果
RT =
0.32050000
勒让德的逼近要求被逼近函数定义在区间【-1,1】上，勒让德多项式也可以通过递推来定义：
它们之间满足如下的正交关系：
在实际应用中，可以根据所需的精度来截取有限项数。勒让德级数中的系数由下式确定：
在MATLAB中编程实现的勒让德逼近法函数为：Legendre。
功能：用勒让德多项式逼近已知函数。
调用格式：f= Legendre（y，k）或f= Legendre（y，k，x0）
分别建立ac，cd，db三段对应的弯矩方程：
M1=Fay*x-q*x.^2/2;0≦x≦L/2
M2=Fby*(L-x)-M0;L/2≦x≦3L/4
M3=Fby*(L-x);3L/4≦x≦L
第三步：建立挠曲轴近似微分方程并积分
建立挠曲轴近似微分方程d2Y/dx2=M(x)/EI
对M/EI积分，得转角A，再做一次积分，得挠度Y，每次积分都有一个待定的积wk.baidu.com常数。
end
R=RT(k+1,k+1);
3）应用romberg数值积分求解
源程序
apply3.m
F=inline('exp((-x.^3)./2)./(sqrt(2*pi))');
[RT,R,wugu,h]=romberg(F,0,1,1.e-8,13)%调用romberg函数计算F的近似积分值symsx;
A=∫(M)*dx/(E*I)+Ca= A0（x）+Ca，此处设A0（x）=cumtrapz(M)*dx/(E*I)
Y=∫(A)*dx+Cy=∫A0(x)*dx+Ca*x+Cy,此处设Y0（x）=cumtrapz(A0)*dx
第四步：确定相应的积分常数
Ca，Cy由边界条件Y(0)=0，Y(L)=0确定
f = 0.5935
从逼近结果上来看，函数的准确值为1/(2-0.5)=0.6667.
④流程图如下：
第一题流程图：
是
否
是
否
是
第二题流程图：
二．编程解决以下科学计算和工程实际问题。
1）简支梁受左半均匀分布载荷q及右边L/4处集中力偶M0作用（如下图1-1），求其弯矩、转角和挠度。设L=2m，q=1000N/m，M0=900N*m，E=200*109N/m2,I=2*10-6m4.
当n=4时，可以查出6阶牛顿-科茨公式的科茨系数 为
, , ,
,
根据牛顿-科茨公式课的六阶牛顿-科茨公式为
2)NewtonCotes函数流程图
否
是
否
是
3)牛顿-科茨数值积分的主程序
NewtonCotes.m
functionI = NewtonCotes(f,a,b,eps)
%牛顿-科茨公式求函数f在区间[a,b]上的定积分
其中，y为已知函数；
k为逼近已知函数所需项数；
x0为逼近点的x坐标：
f为求得的勒让德逼近多项式或在x0处的逼近值。
2源程序：
第一个源程序：
function f=Legendre(y,k,x0)
%用勒让德多项式逼近已知函数
%已知函数：y
%逼近已知函数所需项数：k
%逼近点的x坐标：x0
%求得的勒让德逼近多项式或在x0处的逼近值：f
2）romberg函数流程图
是
否
2)romberg数值积分的主程序
romberg.m
function[RT,R,wugu,h]=romberg(fun,a,b,wucha,m)
%fun为被积函数f(x)
%a为积分下限
%b为积分上限
%m为romberg积分表中行的最大数目
%wucha是两次相邻迭代值得绝对误差限，即|T(j+1,j+1)-T(j+1,j)|<wucha
1）解题思路
这是一个迭代方程组，初始值为x(0)=0,y(0)=0,作为迭代迭代点对(x(i),y(i))的起始点我们可以利用循环迭代将迭代出来的随机点存放在x,y矩阵之中，并调用plot命令将所有点画出来，由于点过于密集，彼此之间的距离很近以致宏观上难以区分，便看似一条连贯的图线。
2）求解主程序
源程序
设函数 闭区间 上有n+1个等距节点 处有定义，以 为节点，对被积函数 建立n次拉格朗日插值多项式。
,则
所以:
其中 ， ，
令 ，则
由此可得:
此公式称为牛顿—科茨求积分公式， 称为求积分的节点， 称为求积分的系数， 称为Cotes系数。
其中梯形公式对应于牛顿-科茨公式n=1的情况，辛普森公式对应于其n=2的情况。且当 时，n越大，积分精度越高：但是n>6时，积分的稳定性得不到保证，因此本题只采用n=6时的牛顿-科茨公式对待解问题进行计算。
Fby=(q*L^2/8+M0)/L;Fay=q*L/2-Fby;%求支反力
图1-1
1解题思路：
首先对简支梁进行受力分析，受力分析图（如下图1-2）所示：
图1-2
第一步：计算支反力
设支座a和b处的支反力分别为Fay和Fby，则据∑Ma=0，∑Fy=0得到平衡方程为：
Fby=(q*L^2/8+M0)/L
Fay=q*L/2-Fby
第二步：建立弯矩方程
以截面c，d为分界面，将梁划分为ac，cd，db三段
Y(0)= Y0(0)+Cy=0
Y(L)= Y0(L)+Ca*L+Cy=0
即[0 1；L 1] [Ca；Cy]=[-Y0（0）；-Y0（L）]
[Ca；Cy]=[0,1;L,1]\[-Y0(0);-Y0（L)];
第五步：根据计算结果绘制弯矩、转角以及挠度图形
②源程序：
L=2;q=1000;M0=900;E=200e9;I=2e-6;%输入已知参数
%函数名
%积分下限a
%积分上限b
%积分I
%精度eps
ifnargin<3
error
return;
else
ifnargin==3
eps = 1.0e-4;%默认精度
end
end
I = ((b-a)/288)*(19*subs(sym(f),findsym(sym(f)),a)+75*subs(sym(f),findsym(sym(f)),(4*a+b)/5)+50*subs(sym(f),findsym(sym(f)),(3*a+2*b)/5)+50*subs(sym(f),findsym(sym(f)),(2*a+3*b)/5)+75*subs(sym(f),findsym(sym(f)),(a+4*b)/5)+19*subs(sym(f),findsym(sym(f)),b));
syms t;
P(1:k+1)=t;
P(1)=1;
P(2)=t;
c(1:k+1)=0.0;
c(1)=int(subs(y,findsym(sym(y)),sym('t'))*P(1),t,-1,1)/2;
c(2)=int(subs(y,findsym(sym(y)),sym('t'))*P(2),t,-1,1)/2;
%RT为romberg积分表
%R为利用romberg求积公式计算被积函数的数值积分值
%wugu为误差估计
%h为最小步长
n=1;h=b-a;wugu=1;x=a;k=0;RT=zeros(4,4);
RT(1,1)=h*(feval(fun,a)+feval(fun,b))/2;
while((wugu>wucha)&(k<m)|(k<4))
k=k+1;h=h/2;s=0；
forj=1:n
x=a+h*(2*j-1);
s=s+feval(fun,x);%计算
end
RT(k+1,1)=RT(k,1)/2+h*s;n=2*n;
fori=1:k
RT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1);
end
wugu=abs(RT(k+1,k)-RT(k+1,k+1));%wugu为误差估计
wR与wugu的绝对误差精确值wR1=3.0131e-010，即两者很接近。
4） 流程图
勒让德逼近编程实现“勒让德逼近”，并用勒让德公式（取6项）逼近函数：1/(2-x)，并求出当x=0.5时的函数值。
1算法说明：
当一个连续函数定义在区间[-1，1]上时，它可以展开成勒让德级数。即：
其中Pn(x)为n次勒让德多项式
f=c(1)+c(2)*t;
for i=3:k+1
P(i)=((2*i-3)*P(i-1)*t-(i-2)*P(i-2))/(i-1);
c(i)=int(subs(y,findsym(sym(y)),t)*P(i),t,-1,1)/2;
f=f+c(i)*P(i);
if(i==k+1)
if(nargin==3)%输入参数个数为三
4)应用六阶牛顿-科茨公式，求定积分
源程序：
apply1.m
symsx;
a=0;b=pi/2;
I=exp(sin(x));
s=NewtonCotes('I',a,b)
运行结果：
s=
1.2337
5)函数应用流程图
6)
7)
2、给出迭代方程
先编写求解方程的函数文件，然后调用该函数文件求30000个点上的x、y,最后在所有的 坐标处标记一个点（不要连线）会出图形，这种图形又称为埃农（Henon）引力线图，它将迭代出来的随机点吸引到一起，最后得出貌似连贯的引力图线。