上面的網(wǎng)址不知道什么時候就打不開了，趕緊保存一份，要不想看都看不到了。

什么是流函數(shù)，什么是位函數(shù)(勢函數(shù))，可以自己搜索。

說說我這里的應(yīng)用場景。

空間放一些電荷，我們能夠算出任意一點的電場強度——一個矢量，現(xiàn)在，我們能不能通過這些矢量來求出電場的等值線，即我可以通過給不同的值設(shè)置不同的顏色，最終得到電場線。這樣得到的電場線不會有偏差，而用矢量一步一步畫出來的電場線必然會有偏差。

Q:?How?to?compute?a?streamfunction?

? Kirill Pankratov, Ph. D. (kirill@plume.mit.edu)

Department of Earth, Atmospheric & Planetary Sciences,

Massachusetts Institute of Technology, Cambridge, MA, 02139

(posted on comp.soft-sys.matlab, 1994-03-07)

Hi,

Although this question appears for the first time?to?my knowledge I believe it can be of rather general interest.

The streamfunction (and also the velocity potential) are rather fun-damental concepts for anybody who is dealing with fluid dynamics or geosciences, like power spectrum for signal processing.

The Mathwolks (I mean the folks from The Mathworks) could be surprised?how?many of the MATLAB users are oceanographers, meteorologists, geophysisists and not only electrical engineers or signal processors (no offence for the latter). So probably The Mathwork should learn at least some fundamentals about it and include the most general operations from these fields in its library.

Does anybody agree with this? ...

Anyway, here I present the program "flowfun" which calculates the streamfunction psi and the velocity potential phi (that is the velocity field is?a?scalar and vector product of the gradient operator and the potential and the streamfunction correspondingly):u?=?d(phi)/dx,???v?=?d(phi)/dy??for?potential?flows,

u?=?-d(psi)/dy,??v?=?d(psi)/dx??for?solenoidal?flows.These potentials are computed by integrating the velocities given by matrices u, v using Simpson rule summation.?To?do this the routine "cumsimp" is added. This routine is actually of much more general application, than the "flowfun" - is is general cumulative integrator -?a?little bit similar?to?"cumsum" but much more accurate (for continuous functions) and starting from zero. I think MATLAB should have had something like this long ago it would be useful for anybody dealing with continuous fields. See HELP for additional information.

Now back?to?"flowfun".

To?get the feeling what it is all about try the following script:e?=?2;?g?=?1;

[x,y]?=?meshgrid(0:20,0:15);??%?This?makes?regular?grid

u?=?e*x-g*y;??????????????????%?Linear?velocity?field

v?=?g*x-e*y;

[phi,psi]?=?flowfun(u,v);??%?Here?comes?the?potential?and?streamfun.

contour(phi,20,'--r')???%?Contours?of?potential

hold?on

contour(psi,20,'-g')????%?Contours?of?streamfunction

quiver(u,v,'w')?????????%?Now?superimpose?the?velocity?field

%?Now?see?the?meaning?of?these?potentials?

If?you?want?the?streamfunction?only,?use

psi?=?flowfun(u,v,'-')

(or?psi?=?flowfun(u,v,'psi')?or?psi?=?flowfun(u,v,'stream')?).I would appreciate any comments and suggestions about these routines. Regards, Kirill.

And now the code itself (flowfun.m and cumsimp.m).==================================??save?as???flowfun.m??=========

function??[phi,psi]?=?flowfun(u,v,flag)

%?FLOWFUN??Computes?the?potential?PHI?and?the?streamfunction?PSI

%?????of?a?2-dimensional?flow?defined?by?the?matrices?of?velocity

%?????components?U?and?V,?so?that

%

%???????????d(PHI)????d(PSI)??????????d(PHI)????d(PSI)

%??????u?=??-----??-??-----?,????v?=??-----??+??-----

%????????????dx????????dy??????????????dx????????dy

%

%?????For?a?potential?(irrotational)?flow??PSI?=?0,?and?the?laplacian

%?????of?PSI?is?equal?to?the?divergence?of?the?velocity?field.

%?????A?non-divergent?flow?can?be?described?by?the?streamfunction

%?????alone,?and?the?laplacian?of?the?streamfunction?is?equal?to%?????vorticity?(curl)?of?the?velocity?field.

%?????The?stepsizes?dx?and?dy?are?assumed?to?equal?unity.

%???[PHI,PSI]?=?FLOWFUN(U,V),?or?in?a?complex?form

%???[PHI,PSI]?=?FLOWFUN(U+iV)

%?????returns?matrices?PHI?and?PSI?of?the?same?sizes?as?U?and?V,

%?????containing?potential?and?streamfunction?given?by?velocity

%?????components?U,?V.

%?????Because?these?potentials?are?defined?up?to?the?integration

%?????constant?their?absolute?values?are?such?that

%?????PHI(1,1)?=?PSI(1,1)?=?0.

%?????If?only?streamfunction?is?needed,?the?flag?can?be?used:

%???PSI?=?FLOWFUN(U,V,FLAG),?where?FLAG?can?be?a?string:

%?????'-',?'psi',?'streamfunction'?(abbreviations?allowed).

%?????For?the?potential?the?FLAG?can?be??'+',?'phi',?'potential'.

%??Uses?command?CUMSIMP?(Simpson?rule?summation).

%??Kirill?K.?Pankratov,?March?7,?1994.

%?Check?input?arguments?.............................................

issu=0;?issv=0;?isflag=0;????%?For?input?checking

isphi?=?1;?ispsi?=?1;????????%?Is?phi?and?psi?to?be?computed

if?nargin==1,?issu?=?isstr(u);?end

if?nargin==2,?issv?=?isstr(v);?end

if?nargin==1&~issu,?v=imag(u);?end

if?issv,?flag?=?v;?v?=?imag(u);?isflag?=?1;?end

if?nargin==0|issu????????????%?Not?enough?input?arguments

disp([10?'??Error:?function?must?have?input?arguments:'...

10?'??matrivces??U?and?V??(or?complex?matrix?W?=?U+iV)'?10?])

return

end

if?any(size(u)~=size(v))?????%?Disparate?sizes

disp([10?'??Error:?matrices?U?and?V?must?be?of?equal?size'?10])

return

end

if?nargin==3,?isflag=1;?end

u?=?real(u);

%?Check?the?flag?string?.?.?.?.?.?.?.?.

Dcn?=?str2mat('+','potential','phi');

Dcn?=?str2mat(Dcn,'-','streamfunction','psi');

if?isflag

lmin?=?min(size(flag,2),size(Dcn,2));

flag?=?flag(1,1:lmin);??A?=?flag(ones(size(Dcn,1),1),1:lmin)==Dcn(:,1:lmin);

if?lmin>1,?coinc?=?sum(A');?else,?coinc?=?A';?end

fnd?=?find(coinc==lmin);

if?fnd~=[],?if?fnd<4,?ispsi=0;?else,?isphi=0;?end,?end

end

phi?=?[];????????%?Create?output

psi?=?[];

lx?=?size(u,2);??%?Size?of?the?velocity?matrices

ly?=?size(u,1);

%?Now?the?main?computations?.........................................

%?Integrate?velocity?fields?to?get?potential?and?streamfunction

%?Use?Simpson?rule?summation?(function?CUMSIMP)

%?Compute?potential?PHI?(potential,?non-rotating?part)

if?isphi

cx?=?cumsimp(u(1,:));??%?Compute?x-integration?constant

cy?=?cumsimp(v(:,1));??%?Compute?y-integration?constant

phi?=?cumsimp(v)+cx(ones(ly,1),:);

phi?=?(phi+cumsimp(u')'+cy(:,ones(1,lx)))/2;

end

%?Compute?streamfunction?PSI?(solenoidal?part)

if?ispsi

cx?=?cumsimp(v(1,:));??%?Compute?x-integration?constant

cy?=?cumsimp(u(:,1));??%?Compute?y-integration?constant

psi?=?-cumsimp(u)+cx(ones(ly,1),:);

psi?=?(psi+cumsimp(v')'-cy(:,ones(1,lx)))/2;

end

%?Rename?output?if?need?only?PSI

if?~isphi&ispsi&nargout==1,?phi?=?psi;?end

===========================?end??flowfun.m?======================

============================?save?as??cumsimp.m?=================

function??f?=?cumsimp(y)

%?F?=?CUMSIMP(Y)????Simpson-rule?column-wise?cumulative?summation.

%???????Numerical?approximation?of?a?function?F(x)?such?that

%???????Y(X)?=?dF/dX.??Each?column?of?the?input?matrix?Y?represents

%???????the?value?of?the?integrand??Y(X)??at?equally?spaced?points

%???????X?=?0,1,...size(Y,1).

%???????The?output?is?a?matrix??F?of?the?same?size?as?Y.

%???????The?first?row?of?F?is?equal?to?zero?and?each?following?row

%???????is?the?approximation?of?the?integral?of?each?column?of?matrix

%???????Y?up?to?the?givem?row.

%???????CUMSIMP?assumes?continuity?of?each?column?of?the?function?Y(X)

%???????and?uses?Simpson?rule?summation.

%???????Similar?to?the?command?F?=?CUMSUM(Y),?exept?for?zero?first

%???????row?and?more?accurate?summation?(under?the?assumption?of

%???????continuous?integrand?Y(X)).

%

%????See?also?CUMSUM,?SUM,?TRAPZ,?QUAD

%??Kirill?K.?Pankratov,?March?7,?1994.

%?3-points?interpolation?coefficients?to?midpoints.

%?Second-order?polynomial?(parabolic)?interpolation?coefficients

%?from??Xbasis?=?[0?1?2]??to??Xint?=?[.5?1.5]

c1?=?3/8;?c2?=?6/8;?c3?=?-1/8;

%?Determine?the?size?of?the?input?and?make?column?if?vector

ist?=?0;?????????%?If?to?be?transposed

lv?=?size(y,1);

if?lv==1,?ist?=?1;?y?=?y(:);?lv?=?length(y);?end

f?=?zeros(size(y));

%?If?only?2?elements?in?columns?-?simple?sum?divided?by?2

if?lv==2

f(2,:)?=?(y(1,:)+y(2))/2;

if?ist,?f?=?f';?end???%?Transpose?output?if?necessary

return

end

%?If?more?than?two?elements?in?columns?-?Simpson?summation

num?=?1:lv-2;

%?Interpolate?values?of?Y?to?all?midpoints

f(num+1,:)?=?c1*y(num,:)+c2*y(num+1,:)+c3*y(num+2,:);

f(num+2,:)?=?f(num+2,:)+c3*y(num,:)+c2*y(num+1,:)+c1*y(num+2,:);

f(2,:)?=?f(2,:)*2;?f(lv,:)?=?f(lv,:)*2;

%?Now?Simpson?(1,4,1)?rule

f(2:lv,:)?=?2*f(2:lv,:)+y(1:lv-1,:)+y(2:lv,:);

f?=?cumsum(f)/6;??%?Cumulative?sum,?6?-?denom.?from?the?Simpson?rule

if?ist,?f?=?f';?end?????%?Transpose?output?if?necessary

=============================?end??cumsimp.m?=================

總結(jié)

以上是生活随笔為你收集整理的matlab计算流函数,hanyeah的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。