Fortran语言编程练习题

平面驻点流动方程数值解

1问题描述

自编程序完成平面驻点流动方程：的数值解。

边界条件为

式中：。

2求解过程

由于上述方程是非线性方程，可采用MATLAB 7.0软件来求解，步骤如下：

第一步：将方程化为一阶常微分方程组。

边界条件为

第二步：建立ode.m和lbc.m两个M文件。

ode.m文件程序如下：

function dfdx=ode(x, f)

dfdx=[f(2);f(3);-f(1)*f(3)+f(2)^2-1];

lbc.m 文件程序如下：

function res=lbc(f0,finf)

res=[f0(1);f0(2);finf(2)-1];

第三步：求解方程。

在MATLAB 7.0工作窗口输入程序：

infinity=4;

solinit=bvpinit(0:0.4:infinity,[0 0 0]);

sol=bvp4c(@ode,@lbc,solinit);

x=0:0.4:infinity

f=deval(sol,x)

plot(x,f(1,:),'ob',x,f(2,:),'rp',x,f(3,:),'b*') /*绘图命令*/

xlabel('轴\it \eta');ylabel('轴\it \phi')

legend('平面驻点流动\phi曲线','平面驻点流动d\phi/d\eta曲线','平面

驻点流动d^2\phi/d\eta^2曲线')

title('平面驻点流动的数值解')

3结果分析

平面驻点流动数值解的计算结果见表1，图1是平面驻点流动数值解的散点图。从表和图中可以看出，从开始呈线性增长，随着的增加，偏离斜直线，当以后渐近于1。在左右，，即此时粘性流动的流速已接近

势流流速，只差百分之一。

00.40000.8000 1.2000 1.6000 2.000 2.4000 2.8000 3.2000

00.08800.31240.62200.9798 1.3620 1.7553 2.1530 2.5523

00.41450.68590.84670.93240.97320.99060.99710.9992

1.23250.84630.52510.29380.14730.06580.02600.00910.0028

表1 平面驻点流动数值解的计算结果

图1 平面驻点流动数值解散点图