Matlab编程以及接口

Matlab编程的初步知识：

1.M脚本文件

Matlab命令集合；

无需输入和输出变量；

与其它脚本文件之间变量透明，即共同使用Matlab的工作空间。

比如

例1：

x = (1:10);

y =sin(x);

plot(x,y)

例2：

A = magic (4);

B = A';

[C,D = xtimesAB(A,B);

2.函数

同样也是完成具体任务的指令集合，但是用一个名字封装起来，变量对外不透明，需要借助输入变量提供数据，输出变量给出结果。执行完毕后，所用的内存全部释放给Matlab。这样的命令集合体称为函数，封装的名字称为函数名，输入变量和输出变量在函数名前后指定。比如例2的xtimesAB(A,B)定义如下：

function [C,D] =xtimesAB(A,B)

%

% This function can tell the products of e-e and V-C

%

C = A.*B;

D = A*B;

end

其中，函数名为xtimesAB，输入变量为A,B；输出变量为C,D。函数名下面%开头的部分为注释内容。在Matlab环境下，可以通过help xtimesAB来显示。

还有一种简单的函数，即inline函数，其特点是随用随定义。比如

1.>> myfun = ‘1+log(r) ‘;

2.>> myfuni=inline(myfun,’r’)

3.>>a=feval(myfuni,10)

4.结果a = 3.3026

我们用得最多的，就是像xtimesAB这样的M函数。与上面这个inline函数对应的M函数为

function y=myfun(r)

y=1+log(r);

该函数结尾不含有end，即Matlab不要求必须有end。

使用时，在Matlab命令环境下，直接书写函数名，并给出输入变量，即可以工作。但需要注意的是，Matlab调用的是文件名而不是函数名。但一般两者取一样的名字。因此，关键是文件名。

M函数的流程控制简介如下：

I. 循环

5.while

a)while condition

b)命令语句

c)end

6.for

a)for i=1:n

b)语句

c)end

ii.条件

1.if

a)if condition 1

b)语句

c)elseif condition2

d)语句

e)else

f)语句

g)end

2.switch

a)switch 表达式

b)case case1

c)语句

d)case case2

e)语句

f)otherwise

g)语句

h)End

iii.跳过：continue

iv.跳出：break

2009.9.29

绘图补充1. EZ系列绘图函数

ezplot('cos(x)')

ezplot('cos(x)', [0, pi])

ezplot('1/y-log(y)+log(-1+y)+x - 1')

ezplot('x^2 - y^2 - 1')

ezplot('x^2 + y^2 - 1',[-1.25,1.25]); axis equal

ezplot('x^3 + y^3 - 5*x*y + 1/5',[-3,3])

ezplot('x^3 + 2*x^2 - 3*x + 5 - y^2')

ezplot('sin(t)','cos(t)')

ezplot('sin(3*t)*cos(t)','sin(3*t)*sin(t)',[0,pi])

ezplot('t*cos(t)','t*sin(t)',[0,4*pi])

f = @(x)cos(x)+2*sin(x);

ezplot(f)

ezplot(@humps)

其它还有

ezcontour, ezcontourf, ezmesh, ezmeshc, ezplot3, ezpolar, ezsurf, ezsurfc Maple

例题1 用二分法计算2

介于 1- 2

先在()1,2内折半，112x =，则2

9

24x =

>； 在11,12?? ???

内折半，114

x =，则2

25

216

x =

<； 在111,142?? ???

内折半，318

x =，得2

121

264

x =

< 在311,182?? ???

内折半，...... 7116x =. (2)

5292256

x => ............

M = 2 a = 1 b = 2 k = 0;

while b-a > eps x = (a + b)/2; if x^2 > M b = x else a = x end k = k + 1; end

共计52步，得到 1.4142。 慢但很可靠。

上面问题也可以写成

()22f x x =-

因此，上面的计算过程就是求根过程：()0f x =。可以用牛顿法进行计算，即

()()

1'n n n n f x x x f x +=-

即

212

2n n n n

x x x x +-=-

一般地，计算M 的算法，

21212n n n n

n n x M

x x x M x x +-=-

??=+ ???

M=2; x=1; xprev=2;

while abs(x - xprev) > eps*abs(x) xprev = x; x = 0.5*(x + M/x); end

6步迭代即计算完毕，比上面快，但有条件限制

1）连续函数； 2）导数方便可求； 3）初值距离解要很近。

牛顿法的特例： 对于

当1n x 和n x 围绕a 振荡不收敛时，有

此时的函数满足

则

当 a=2, 则ezplot('sign(x-2)*sqrt(abs(x-2))',0,4)

导数就改点切线的斜率，用割线斜率代替切线斜率，称为割线法，即

a = 1

b = 2

k = 0;

while abs(b-a) > eps*abs(b)

c = b - (b - a)/(1-f(a)/f(b));

k = k + 1;

a = b;

b = c;

end

function y=f(x)

y = x^2 - 2;

这里，f(x)就是M函数。在M脚本文件中，不能直接使用这个定义，需要在单独文件中进行定义。

该函数也可以用Fortran编写，然后在Matlab内调用，即Fortran得mex文件。

mex文件

dll

Matlab可以和C/Fortran进行混合编程。

混合编程有两个方面：C/Fortran调用Matlab库；二是，Matlab内调用C/Fortran的函数。

mex文件，就是把C/Fortarn程序，编译成Matlab动态库文件，然后在matlab内调用。这样Fortran程序可以放在matlab中使用，对于已有Fortran软件来说，是一件方便之举。

Fortran程序书写完毕，存入文件后，就可以进行编译。编译方法就是：

mex mycode.f

之前，计算机内要装fortran编译器。

Fortran程序中，自己的函数一定要用subroutine来定义；此外，还要加入一个matlab的接口函数，

subroutine mexFunction(nlhs, plhs, nrhs, prhs)

程序f.f清单如下

subroutine f(x,y)

real*8 x,y

y = x*x - 2

end

C The gateway routine

subroutine mexFunction(nlhs, plhs, nrhs, prhs)

integer mxGetM, mxGetN, mxGetPr

integer mxIsNumeric, mxCreateDoubleMatrix

integer plhs(*), prhs(*)

integer*4 x_pr, y_pr

integer nlhs, nrhs

integer m,n,size

C Check for proper number of arguments.

if(nrhs .ne. 1) then

call mexErrMsgTxt('USAGE: y=f(x).')

endif

C Get the size of the input array.

m = mxGetM(prhs(1))

n = mxGetN(prhs(1))

size = row*col

x_pr = mxGetPr(prhs(1))

plhs(1) = mxCreateDoubleMatrix(m,n,0)

y_pr = mxGetPr(plhs(1))

C Call the computational subroutine.

call f(%val(x_pr),%val(y_pr))

return

end

编译命令：

>> mex f.f

用Matlab内的文本编辑器书写上面程序。注意，matlab兼容fortran77语言规范，即（1）每行最多含72个字符，而且前6个位置要空出来，每行程序语句起始于第7列；

（2）当某一行为注释文字，就在第一列键入* 或者C;

（3）2-5用于行号标注；

（4）当公式太长，一行容不下，则折转到从下一行，并从第7个字符位置开始。同时，在第6个字符位置，键入任何字符，表示这是续上一行。

例题2 封闭空间辐射换热的有效辐射法

1．基本概念

（1）漫表面与灰表面

发射率和吸收率不随发射方向改变的表面，称为漫射表面；发射率和吸收率不随波长改变的表面，称为灰表面；兼有两者性质的表面，称为灰漫表面。 在半球方向上，灰漫表面辐射强度不变，总辐射功率记为E 。

（2）辐照量( irradiation )

来自其它表面的所有投射辐射能（包括其它表面反射过来的能量）G 。

（3）有效辐射( radiosity )

离开灰漫表面的所有辐射，包括发射和反射两个部分，应为J ；

（4）表面总有效辐射计算

b J E G ερ=+

,ερ表面发射率和反射率。

2．封闭空间内单个表面的辐射

房间是封闭空间的，房间内的空气为透明气体。假定每个墙体表面是等温的灰漫表面后，应用有效辐射概念可以直接进行墙体表面之间的辐射换热计算。这里，离开表面的有效辐射和到达表面的有效投射，能流强度分布均匀。

到达表面的总热量记为i q ，在热平衡状态下，应等于离开表面i 的辐射量，即

()i i i i q A J G =-

i J 为离开表面的辐射能和反射能，即有效辐射（radiosity ）为，

i i i i J E G ρ=+

i G 为到达表面的辐射能，即有效辐照（irradiation ）。

对于不透明的灰漫表面（ opaque, diffuse, gray surface ）,有11i i i ραε=-=-，则有效辐射为

()1i i bi i i J E G εε=+-

则

1i i bi

i i i i J E q A J εε??

-=- ?-?

?

或

()1bi i

i i i i

E J q A εε-=

-

此公式描述了离开表面的净辐射量，与表面的黑体辐射与有效辐射之差有关系，可直观地表示电阻为()1i i i A εε-的网络单元，其驱动势差为黑体辐射功率与有效辐射功率之差。

3． 两个表面之间的辐射换热

已知两个灰漫表面的有效辐射，就可以采用黑体辐射换热计算方法确定这两个表面之间的辐射换热，因此i J 是关键的中间量。为确定此值，需要同时计算所有表面之间的辐射换热计算等式。如同黑体之间辐射换热计算一样，计算之前也需要确定表面之间的角系数ij F ，即表面i 到j 的角系数。

封闭空间内，角系数具有以下性质 （1） 互易关系：

ij j ji i F A F A =

（2） 累加原则：

1

1N

ij

i F

==∑

表面j 到表面i 辐射量为

'ij ij j j q F A J =

表面i 到表面j 辐射量为

'ji ji i i q F A J =

表面j 到表面i 的净辐射量为

()

1

j i

ij j ij

J J q A F --=

表示为电阻为()

1

j ij

A F -的网络单元, 驱动力为()

j i J J -。

4． 封闭空间内所有表面之间的辐射换热

封闭空间有N 个表面，根据上式，表面j 到其余表面的净辐射量为

()

1

1

1

N N

j i

j ij i i j ij

J J q q A F -==-==∑∑

因此， 对于封闭空间，存在一线性方程组，

[]A ?=J B

可以计算所有表面的有效辐射[]12,,,N J J J =J ，[]12,,,N B B B =B 。

5．表面的类型

1）温度已知表面，则4

bi i E T σ=，计算其热流量i q

()1bi i

i i i i

E J q A εε-=

-

2）热流i q 已知的表面，则4

bi i E T σ=联合下式，计算表面温度i T 。

()1bi i i i i i E J q A εε=+-

3）特殊情况

A ）表面为绝热表面（reradiating surface ），则0i q =；

B ）黑体表面，则bi i E J =

6．表面有效辐射的计算

（1）表面温度已知

()()1

1N

bi i

i ji i j j i i i E J A F J J A εε=-=--∑

()()1

1N

bi i i ji i j i j E F J J J εε==-?-+∑

()

111N N bi i i ji i ji j i j j E F J F J J εε==??

=-?-+????∑∑

()

11N bi i i i ji j i j E J F J J εε=??

=-?-+????

∑

()()1

111N

bi i i i i i ji j j E J F J εεεε==+---????∑

1

N

i ii i ji j j b a J a J ==+∑

11i

ii i

a εε-=+

，1i

ij ji i

a F εε-=-

，4i bi b E T σ==

（2） 表面热流已知

()1

N

i i ji j i j q A F J J ==-∑

()1

N

i

i ji j j i q J F J A ==-∑ 1ii a =，ij ji a F =-，i

i i

q b A =

（3） 特殊情况

i. 黑体表面：1i ε=，计算格式采用（1）； ii.

绝热表面：0i q =，计算格式采用（2）。

7．编制Matlab 程序

a) 函数名及其调用方式：

J=radiosity( [type1, T1 or q1, eps1, Area1, F12, F13, …],

[type2, T2 or q2, eps2, Area2, F23, …], [type3, T3 or q3, eps3, Area3, …], [……])

其中，每一个参数代表一面墙体，每一面墙体用下面参数描述，并用[]括在一起： 墙体类型：type

1 等温表面

2 等热流表面

对应的物理参数为：T 温度

Q 热流

表面辐射率（eps ）

i ε

面积（Area ）

i A

对其它房间的角系数F ，

比如三个表面的角系数矩阵为 F11 F12 F13 F21 F22 F23

F31 F32 F33

输入右上三角

F12 F13

F23

b) 单元变量

单元（cell ）是Matlab 的一个数据类型，其内部含有若干元素，元素本身可以是任何数据类型，数，矩阵，字符串等。其初始化的方法如下，

S = {

[type1, T1 or q1, eps1, Area1, F12, F13]; [type2, T2 or q2, eps2, Area2, F23]; [type3, T3 or q3, eps3, Area3]

}

因此， 生成向量()3,1S ，其长度为()length S ＝3，每个元素也为向量。

1S ＝[type1, T1 or q1, eps1, Area1, F12, F13]； 2S ＝[type2, T2 or q2, eps2, Area2, F23]； 3S ＝[type2, T2 or q3, eps3, Area3]；

c) 变参数函数

Varargin 和varargout 是Matlab 的一对系统哑元, 专门用于传递数目变化的数据，其使用方法为

function varargout = radiosity (n)

我们只想了解varargin ，即

function y = radiosity (varargin)

当使用radiosity 函数进行计算时，不同题目，墙体数目会发生变化。具体的参数数目可以从系统变量nargin 中读取。参数的数目，是指逗号之间的参数个数，比如

J=radiosity( [type1, T1 or q1, eps1, Area1, F12, F13],

[type2, T2 or q2, eps2, Area2, F23], [type3, T3 or q3, eps3, Area3] ）

参数为三个，即 nargin 内容为3。

也可以用length(varargin)返回参数数目。

每个墙体的参数包含在varargin 每一个参数的元素内，即 第一个墙体表面为 varargin{1}

表面类型为varargin{1}（1）

表面温度为varargin{1}（2）

表面发射率varargin{1}（3）

表面面积varargin{1}（4）

角系数varargin{1}（5：end）d) 程序

function J=radiosity(varargin)

sigma=5.67e-8;

m=nargin;

eps=zeros(m,1);

J=zeros(m,1);

type=zeros(m,1);

A=zeros(m,1);

F=zeros(m,m);

a=zeros(m);

b=zeros(m,1);

for i=1:m

type(i)=varargin{i}(1);

phys(i)=varargin{i}(2);

eps(i)=varargin{i}(3);

A(i)=varargin{i}(4);

if(i

F(i,i+1:m)=varargin{i}(5:end);

end

end

for i=1:m

f=(1-eps(i))/eps(i);

switch type(i)

case 1

a(i,i)=1+f;

T=phys(i);

b(i)=sigma*T^4;

case 2

a(i,i)=1;

q=phys(i);

b(i)=q/A(i);

end

if(i

F(i+1:m,i)=F(i,i+1:m)'./A(i+1:m)*A(i);

end

for j=1:m

if i==j

continue

end

switch type(i)

case 1

a(i,j)=-f*F(i,j);

case 2

a(i,j)=-F(i,j);

end

end

F

a

b

J = a\b end

（1）例题13.6（不含黑体表面） （2）例题13.4（含黑体表面）

三面：

T1=1000; ε1=0.9; A1=10; F12=0.39; F13=0.61; T2=600; ε2=0.5; A2=15; F23=0.41； T3=300; ε3=1; A3=20;

计算： J

（4） 习题13.86不含黑体

长方体，10 X 6 X 4。

天棚（10 x 6）ε1=0.8, T1=40 ； 地板（10 x 6）ε2=0.9, T2=50 ； 前立面（10 x 4）ε3=0.7, T3=15 ；

其余面为绝热表面。 计算表面的J 。

10．表面温度或者热流的计算

已知表面温度i T ，则

()4i i i q A T J σ=-

已知表面热流i q ，则

1

i i i i

T q J A =

+

根据上述程序所建立的计算程序如下

作业：编写计算

3的程序。要求

（1）用M 脚本文件； （2）用M 函数； （3）用mex 文件。