Matlab第二、三次上机作业

%文件名：GLAppro.m function[poly,yy,delta]=GLAppro(f,n,a,b,xx) %功能：利用 Gauss Legendre 多项式求函数的最佳平方逼近 %输入：f——被逼近函数；a,b——逼近区间；xx——欲求的逼近点 %n——逼近的 L 多项式的次数（标量时为最高次数，向量时为其所选择的的逼近次数） %输出：poly——所求的逼近多项式系数（降序） ；yy——逼近店的值；delta——逼近误差 N=numel(n); if N>1 id=n+1; else N=n+1;
三：数值试验要求： 1. 试验函数： f ( x ) x cos x , x [0, 4] ；也可自选其它的试验函数。 2. 用所编程序直接进行计算，检测程序的正确性，并理解算法。 3. 分别求二次、三次、 。 。 。最佳平方逼近函数 s( x） 。 4. 分别作出逼近函数 s( x） 和被逼近函数
分别求两次、三次和四次最佳平方逼近函数 s(x) 输入： homework2 输出： 所求的最佳评分那个逼近多项式为（n=2): - 0.1599*x^2 - 0.572*x + 0.8268 误差为:0.45032 所求的最佳评分那个逼近多项式为（n=3): 0.3818*x^3 - 2.45*x^2 + 3.093*x - 0.3948 误差为:0.02393 所求的最佳评分那个逼近多项式为（n=4): 0.08539*x^4 - 0.3014*x^3 - 0.6939*x^2 + 1.531*x - 0.08255 误差为:0.002258
(4)为光滑管摩擦系数计算公式 (5)为粗糙管摩擦系数计算公式 设定λ 初值，带入进行试差，直到计算出的λ 与上一个λ 相差在误差范围内时，试 差结束。 1、算法程序编写 以下为相应算法：
M 文件 1： function [d,u,Re,s,i]=ghshicha(hf,l,v,p,cp) %函数为试差法计算光滑管管径及流速 %输入：hf：压头损失/m l：管长/m v：流量/(m3/h) p：液体密度/(kg/m3) 体粘度/(Pa·s) %输出：d:管径/m u：流速/(m/s) Re：雷诺数 s：摩擦系数 i：迭代次数 s=0.001; for i=1:1000 ss=s; d=(s*v*v*l/(3600*3600*15.41*hf))^(1/5); u=v/(3600*0.78539815*d*d); Re=d*u*p/cp; s=0.0056+0.500/Re^0.32; if s-ss<0.0001 break end end
d= u = ������ Re =
������ ������������������ ������ × ������ 2 4
5
������ ×2������ × ������������ 4
������ × ������ × ������
(1)
(2) (3) (4) (5)
������ = 0.0056 + 0.500/������������ 0.32 ������ = 0.100 × (������ /������ + 68/������������)0.23
wenku.baidu.com
f ( x ) 的曲线图进行比较。
（分别用绘图函数 plot( x 0 ,s( x 0 ))和 fplot(‘xcosx’,[x1 x2,y1,y2])）
解题思路：参照应用数值分析书 P259“算法 7-1” ，利用 Legendre 多项式对 f(x)∈C(a,b)的最 佳平方逼近写出以下算法： M 文件 1：
第二次上机作业
一. 任务： 用 MATLAB 语言编写连续函数最佳平方逼近的算法程序（函数式 M 文件） 。并用此程序进 行数值试验，写出计算实习报告。 二. 程序功能要求： 1. 用 Lengendre 多项式做基，并适合于构造任意次数的最佳平方逼近多项式。
P0 ( x ) 1, P1 ( x ) x 可利用递推关系 Pn ( x ) (2n 1) xPn 1 ( x ) ( n 1) Pn 2 ( x ) / n n 2, 3, .....
2. 程序输入： 3. 程序输出： （1）待求的被逼近函数值的数据点 x 0 (可以是一个数值或向量) （2）区间端点：a,b。 （1）拟合系数： c0 , c1 , c2 ,..., cn （2）待求的被逼近函数值
s( x0 ) c0 P0 ( x0 ) c1 P1 ( x0 ) c2 P2 ( x0 ) cn Pn ( x0 )
乙醇气液平衡拟合
4 3 2 1 0 0.0022 0.0027 乙醇气液平衡拟合 0.0032 0.0037 0.0042
y = -2135.x + 8.058 R² = 0.999
线性 (乙醇气液平衡拟合)
A=8.0582-0.201124=7.857076
B=-2135.3/273=7.8216
调整图象并标注：
1 真值 逼近次数：2 逼近次数：3 逼近次数：4
0
-1
-2
y
-3
-4
-5
0
0.5
1
1.5
2 x
2.5
3
3.5
4
第三次上机作业
一. 任务一： 用 MATLAB 语言编写化工原理流体力学中 d-u-Re-λ 试差法的算法程序（函数式 M 文 件） 。并用此程序进行数值试验，写出计算实习报告。 任务二： 用 MATLAB 语言编写物理化学气液平衡实验数据拟合的算法程序（函数式 M 文件） 。 并用此程序进行数值试验，写出计算报告。 二. 程序功能： 通过编写迭代的算法程序， 使得常规计算过程的人工手算方法得以程序化， 节省时间且 计算准确。 计算示例 1： 一管路总长 70m，要求输水量 30m3/h，输送过程的允许压头损失为 4.5m 水柱。求管径 d。已知水的密度为 1000kg/m3,粘度为 0.001Pa· s； （1）管为光滑管； （2）管内部的绝对粗糙度为 0.2mm。 解： 在此过程中需要用到以下公式：
而利用 MATLAB 编程的计算算法如下： M 文件 1： %文件名：DiscreteOrthoFit.m function[p,pm,pp,delta,yy]=DiscreteOrthoFit(x,y,m,w,xx) %功能：一维离散数据的正交多项式（Schmidt 法）最小二乘拟合 %输入：x,y-初始拟合数据（列向量） ；m-拟合多项式的次数 %w-权系数（可选，默认是 1） ；xx-待求点 %输出：p-m 次多项式的降序系数；pp-对应拟合曲线由正交多项式 Pk（x）组合系数 %pm-其第 k（1<=k<=m+1）行返回对应的 k-1 次正交多项式 Pk（x）的降序系数 %delta-拟合平方误差；yy-待求点数据 if nargin<4||isempty(w);w=ones(length(x),1);else w=w(:);end x=x(:);y=y(:);delta=0; pm=fliplr(eye(m+1,m+1));p=zeros(1,m+1); pm(2,m+1)=-sum(w.*x)/length(x); for k=2:m t1=polyval(pm(k,:),x);t2=polyval(pm(k-1,:),x); a=(x.*w.*t1)'*t1/((t1.*w)'*t1);b=((t1.*w)'*t1)/((t2.*w)'*t2); pm(k+1,:)=conv([1 -a],pm(k,2:m+1))-b*pm(k-1,:);
M 文件 4：
function p=LegPoly(n) %功能：递归法求 n 次 Gauss Legendre 多项式 %输入：n——多项式次数 %输出：p——降幂排列的多项式系数 p0=1;p1=[1,0]; if n==0 p=p0; elseif n==1 p=p1; else p=((2*n-1)*[LegPoly(n-1),0]-(n-1)*[0,0,LegPoly(n-2)])/n; end
M 文件 5：
%文件名：homework2.m f=@(x)(x.*cos(x));n=cell(3,1);n{1}=2;n{2}=3;n{3}=4; a=0;b=4;x0=linspace(a,b)';color=['k','g','b']; y=f(x0);plot(x0,y,'r-','linewidth',1.5);hold on;syms t x; for i=1:3 [poly,py,delta]=GLAppro(f,n{i},a,b,x0);pt=vpa(poly2sym(poly,t),4); poly=simple(subs(pt,t,(2*x-a-b)/(b-a)));poly=vpa(poly,4); disp(['所求的最佳评分那个逼近多项式为（n=' int2str(n{i}) '):']);disp(poly); disp(['误差为:' num2str(delta)]);plot(x0,py,color(i),'linewidth',1.5) end
M 文件 2： function [d,u,Re,s,i]=ccshicha(hf,l,v,p,cp,e)
cp：液
%函数为试差法计算粗糙管管径及流速 %输入：hf：压头损失/m l：管长/m v：流量/(m3/h) p：液体密度/(kg/m3) cp：液体粘 度/(Pa·s) e：粗糙度/mm %输出：d:管径/m u：流速/(m/s) Re：雷诺数 s：摩擦系数 i：迭代次数 s=0.001; for i=1:1000 ss=s; d=(s*v*v*l/(3600*3600*15.41*hf))^(1/5); u=v/(3600*0.78539815*d*d); Re=d*u*p/cp; s=0.100*(e/1000*d+68/Re)^0.23; if s-ss<0.0001 break end end 三、结果检测： 1、 输入：[d,u,Re,s,i]=ghshicha(4.5,70,30,1000,0.001) 输出：d =0.0648u =2.5299Re =1.6384e+005s =0.0163i =3 2、 输入：[d,u,Re,s,i]=ccshicha(4.5,70,30,1000,0.001,0.2) 输出：d =0.0651u =2.5010Re =1.6290e+005s =0.0168i =3 结果显示， 粗糙管由于存在摩擦消耗， 相应的管径要加粗一些才能使输送压头损失保持 在 4.5m。通过编写该算法，日后遇到相同或相似问题可直接调用程序进行计算，不必再人 工手算，提高时效性。
M 文件 2：
function y=myfun(t,f,a,b) %功能：GLAppro 子函数，变换到区间[-1,1] x=(b-a)*t/2+(b+a)/2;y=f(x).*f(x);
M 文件 3：
function f=fLegPoly(t,f1,n,a,b) %功能：GLAppro 子函数，求变换后的积分函数 p=LegPoly(n); x=(b-a)*t/2+(b+a)/2; f=f1(x).*polyval(p,t);
计算示例 2：
方程式转化为：lgP = A + lgP0 −
������������0 ������
命 lgP = y， T = x 构造函数 y = A + lgP0 − ������������0 ������
1
首先用 exel 进行数据预处理，并可得到线性拟合结果，右图：
T/℃ 0 20 40 60 80 100 120 140 160
P/kPa 1.589 5.87 17.89 46.81 108.2 223.6 423.9 747.4 1240
1/T(1/K) 0.003663 0.003413 0.003195 0.003003 0.002833 0.002681 0.002545 0.002421 0.002309
lgP 0.201124 0.768638 1.25261 1.670339 2.034227 2.349472 2.627263 2.873553 3.093422
id=1:N; end delta=quad(@myfun,-1,1,[],[],f,a,b);c=zeros(1,N);poly=zeros(1,id(N)); for k=1:N c(k)=(2*id(k)-1)*quad(@fLegPoly,-1,1,[],[],f,id(k)-1,a,b)/2; delta=delta-c(k)^2*2/(2*id(k)-1); end if nargin==5 t0=(2*xx-a-b)/(b-a);yy=zeros(size(xx)); for k=1:N p=LegPoly(id(k)-1);yy=yy+c(k)*polyval(p,t0); poly(id(N)-id(k)+1:(id(N)))=poly(id(N)-id(k)+1:(id(N)))+c(k)*p; end else for k=1:N p=LegPoly(k-1);poly(N-k+1:N)=poly(N-k+1:N)+c(k)*p; end end