数值计算实例

数值计算

插值

假设需要得到x 坐标每改变0.1 时的y 坐标, 用三次插值方法对机翼断面下缘轮廓线上的部分数据加细, 并作出插值函数的图形.

程序: clear, close all

x=[0,3,5,7,9,11,12,13,14,15]; y=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];

plot(x,y);

xi=0:0.1:15;

yi_cubic=interp1(x,y,xi,'cubic');

plot(x,y,'ro',xi,yi_cubic);

pp=csape(x,y,'second'); v=ppval(pp,xi);

v;

T=(ppval(pp,0.1)-ppval(pp,0))/0.1;

angle=atan(T)*180/pi;

s=v(130:151);

ss=min(s);

图形:

最小二乘拟合

已知空气温度与动力粘度关系如下，进行最小二乘拟合

0℃170.8×10^-4mPa.s

40℃190.4×10^-4mPa.s

74 ℃210.2×10^-4mPa.s

229 ℃263.8×10^-4mPa.s

334℃312.3×10^-4mPa.s

409℃341.3×10^-4mPa.s

481℃358.3×10^-4mPa.s

565℃375.0×10^-4mPa.s

638℃401.4×10^-4mPa.s

750 ℃426.3×10^-4mPa.s

810 ℃441.9×10^-4mPa.s

程序：

>> x=[0 40 74 229 334 409 481 565 638 750 810];

>> y=[170.8 190.4 210.2 263.8 312.3 341.3 358.3 375.0 401.4 426.3 441.9]; >> p=polyfit(x,y,2)

p =

-0.0002 0.4652 172.5460

>> xi=[0:2:810];

>> yi=polyval(p,xi);

>> plot(x,y,'ko-',xi,yi,'k--')

解线性方程组的直接法

某同学做了一个关于气体的热力性质的实验，研究气体在不同的温度和压力情况下熵的变化

S=Aln(T2/T1)–Bln(P2/P1)+C

数学建模：

将实验有关数据代入关系式，得

0.0255 = Aln(600/300)-Bln(1/0.1)+C

0.8496 = Aln(1440/720)-Bln(0.1/0.2)+C

-0.5076 = Aln(360/720)-Bln(0.1/0.2)+C

整理得，

0.69A-2.3B+C=0.0355

0.69A+0.69B+C=0.8596

0.69A-0.69B-C=0.4976

将此方程组看成x的一个列向量进行求解。

程序：

function[x,Aug]=Gaosixiaoqu (A ,b)

n=length(b);

x=zeros(n,1);

c=zeros(1,n+1);

Aug=[A b];

for k=1:n-1

if Aug(k,k)==0

'“矩阵奇异，无解”'

break

end

for i=k+1:n

m=Aug(i,k)/Aug(k,k);

Aug(i,k:n+1)=Aug(i,k:n+1)-m*Aug(k,k:n+1);

end

end

%回代

x(n)=Aug(n,n+1)/Aug(n,n);

for i=n-1:-1:1

x(i)=(Aug(i,n+1)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);

End

在MA TLAB 的命令窗口中输入：

>> A=[0.69 -2.3 1;0.69 0.69 1;0.69 -0.69 -1];

>> b=[0.0355 0.8596 0.4976];

>> [x,Aug]=Gaosixiaoqu(A,b)

x =

0.9835

0.2756

-0.0092

公式应为：

S=0.9835ln(T2/T1)–0.2756ln(P2/P1)-0.0092

牛顿迭代法

问题：

简单蒸馏时，某时刻釜残液量与低沸点组成x 之间的关系式为

⎪⎪⎭

⎫ ⎝⎛--+-=

0011ln ln 11ln

x x x x F

F αα 对于苯-甲苯物系，相当挥发度α=2.5，开始时物系中含苯60%，甲苯40%。试求当蒸馏至残液量为原加料量的一半时残液中的苯含量。 将方程变为： 011ln ln 11ln

f(x)000

=⎪⎪⎭⎫ ⎝⎛--+--

=x x x x F

F αα 将参数代入方程得：06.011ln 5.26.0ln 15.215.0ln

f(x)0

=⎪⎭

⎫

⎝⎛--+--

=x x F F

化简得： 07402.0)l n (x -1ln 5.2f(x)=+-=x ）（

程序：

function[x,k,y]=newton(f,df,x0,tol,N) %Newton 迭代

%f 和df 分别表示函数f(x)及其导函数

%x0表示初值，tol 表示误差限，N 表示最大迭代次数

%输出参数x 表示迭代结果，k 表示迭代次数，y 表示在x 处的函数值f （x ） fprintf('x(%2d)=%10.8f\n',0,x0) for k =1:N

x=x0-f(x0)/df(x0);

err=abs(x-x0); x0=x; y=f(x);

fprintf('x(%2d)=%10.8f\n',k,x) if(err

在命令窗口中

f=inline('2.5*log(1-x)-log(x)+0.7402');