How to compute a streamfunction?
来源:互联网 发布:mysql压缩包下载地址 编辑:程序博客网 时间:2024/06/05 16:15
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 checkingisphi = 1; ispsi = 1; % Is phi and psi to be computedif nargin==1, issu = isstr(u); endif nargin==2, issv = isstr(v); endif nargin==1&~issu, v=imag(u); endif 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 ]) returnendif any(size(u)~=size(v)) % Disparate sizes disp([10 ' Error: matrices U and V must be of equal size' 10]) returnendif nargin==3, isflag=1; endu = 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, endendphi = []; % Create outputpsi = [];lx = size(u,2); % Size of the velocity matricesly = 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 PSIif ~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 vectorist = 0; % If to be transposedlv = size(y,1);if lv==1, ist = 1; y = y(:); lv = length(y); endf = zeros(size(y)); % If only 2 elements in columns - simple sum divided by 2if lv==2 f(2,:) = (y(1,:)+y(2))/2; if ist, f = f'; end % Transpose output if necessary returnend % If more than two elements in columns - Simpson summationnum = 1:lv-2; % Interpolate values of Y to all midpointsf(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) rulef(2:lv,:) = 2*f(2:lv,:)+y(1:lv-1,:)+y(2:lv,:);f = cumsum(f)/6; % Cumulative sum, 6 - denom. from the Simpson ruleif ist, f = f'; end % Transpose output if necessaryMATLAB画流函数
0 0
- How to compute a streamfunction?
- 正在翻译HOW TO THINK LIKE A COMPUTE SCIENTIST ——LEARN PYTHON
- Kgraph-how to compute index file size
- How to Prepare the Windows Azure Compute Emulator
- How to compute the square root of an integer?
- Write a method to compute all permutations of a string
- A program to compute word length of a machine
- the key specified to compute a hash value is expired
- Simple tutorial for using TensorFlow to compute a linear regression
- Manifold SLIC A Fast Method to Compute Content-Sensitive Superpixels
- Compute how many bits to require convert one integer to another
- Compute how many bits to require convert one integer to another
- How to write a Makefile
- How to write a Makefile
- how to write a makefile
- how to become a hacker
- How to Become a Hacker
- How to Read a Paper
- js中两种定时器,setTimeout和setInterval的区别
- zipline的bundle相关数据结构
- ajax 设置Access-Control-Allow-Origin实现跨域访问
- 跨域访问cookie之CORS的完美解决方案
- LeetCode 4. Median of Two Sorted Arrays
- How to compute a streamfunction?
- 机器学习学习笔记(3)----模型评估与选择
- 国足出线的概率
- C# 图像无损压缩
- Linux文件共享原理
- 每天一个Linux命令(20):find命令之exec
- android 点击监听,替换图片
- oracle按照日期自动进行表分区
- SpringMVC有用归纳