function err = errfun(x,D2,pfo,Re,D,sTR2,cotT,R) % ERRFUN - Navier-Stokes error function. % % Syntax: err = errfun(x,D2,pfo) % % This function returns the error in the Navier-Stokes equations, so % to solve the N-S equations means to find x such that errfun(x,...) = 0. % Recover grid size [M,N] = size(D2); % Turn x from a vector into a matrix x = reshape(x,2*M,N); % Pick out the x1=psi and the x2=D2(psi) components x1 = x(1:M,:); x2 = x(M+(1:M),:); % The equation x2 = D2(x1) f1 = D2 * x1 - x2; % Boundary conditions: psi = 0 at inner boudnary and symmetry lines, free flow at outer. f1(1,:) = x1(1,:); f1(:,1) = x1(:,1); f1(:,end) = x1(:,end); f1(end,:) = x1(end,:) - pfo; % The equation D2(D2(psi)) = 0 (only works for Re = 0 at this point) A = D2 - Re./sTR2*((D.dth*x1)*D.dr - (D.dr*x1) * D.dth + ( 2* cotT.* (D.dr*x1) - 2*(D.dth*x1)./R) * eye(D.dr)); f2 = A * x2; % Boundary conditions: D2(psi) = 0 at outer boundary and symmetry lines; at inner % boundary we have D.dr(psi) = 0 f2(1,:) = x1(1,:) - x1(2,:); f2(:,1) = x2(:,1); f2(:,end) = x2(:,end); f2(end,:) = x2(end,:); % Stack up the two error components and return as a vector. err = [f1;f2]; err=err(:); return