使用Chebfun求解Blasius方程(一)

Chebfun的特点:

1. 基于Chebyshev展开，展开项数由机器精度自适应控制;

2. 将符号计算和数值计算结合，以处理数值的速度处理函数；

3. 在Matlab上实现，将Matlab处理向量和矩阵的命令重载，以处理函数和算子；

4. 基于Newton迭代法求解非线性微分方程;

5. 使用自动微分技术计算Frechet导数；

6. Chebop的实现利用了谱方法和惰性求值的思想

7. 能表示具有可去奇点的函数

 

Chebfun仍然在持续开发中，后期对较复杂的求解过程进行了封装，使得用户将更多的精力放在自己的问题上。下面从三个层面使用chebfun系统求解Blasius方程，“麻雀虽小五脏俱全”，以期对Chebfun求解非线性边值问题有较深入的理解。

 

Blasius方程模拟了一个半无限长平板上的二维粘性层流流动，控制方程是

边界条件是

第一个层面求解：

% SolveBlasiusCase1.m
clc;
clear;

infty=100;
N=chebop(0,infty);
N.op=@(x,f) diff(f,3)+f.*diff(f,2)/2;
N.lbc=@(f) [f,diff(f)];
N.rbc=@(f) diff(f)-1;
x=chebfun('x',[0,infty]);
lambda=4;
N.init = x+(1-exp(-lambda*x))/lambda;  
f=N\0;
df=diff(f);
d2f=diff(f,2);

fprintf('Kowarth reports f''''(0) = 0.332057.\n'); 
fprintf('Value computed here is f''''(0) = %7.10f.\n',d2f.vals(1));  

第二个层面求解(将"\"表示的非线性求解方法，即Newton迭代法显式地写出)

% SolveBlasiusCase2.m
clc;
clear;
infty=10;
x=chebfun('x',[0 infty]);
N=chebop(0,infty);
N.op=@(u) diff(u,3)+u.*diff(u,2)./2;
N.lbc=@(u) [u,diff(u)];
N.rbc=@(u) diff(u)-1;

lambda=4;
u = x+(1-exp(-lambda*x))/lambda;  
nrmdu = Inf;
while nrmdu > 1e-10
   J=diff(N,u);
   delta = J\(-N(u));
%    delta =- J\(N*u);
   u = u + delta; 
   nrmdu = norm(delta)
end
du=diff(u);
d2u=diff(u,2);

fprintf('Kowarth reports f''''(0) = 0.332057.\n'); 
fprintf('Value computed here is f''''(0) = %7.10f.\n',d2u.vals(1))

第三个层面求解(将Frechet导数显式地写出)

% SolveBlasiusCase3.m
clc;
clear;
infty=10;
N=@(u) diff(u,3)+u.*diff(u,2)./2; % ODE part of the nonlinear BVP
A=chebop(0,infty);                % Operator on the interval [0,infty]
x=chebfun('x',[0 infty]);         % The chebfun "x" on the problem interval
lambda=4;
u = x+(1-exp(-lambda*x))/lambda;  % initial guess
nrmv = 1; 
while nrmv > 1e-10                % Newton iterations
   A.op=@(v) diff(v,3)+u.*diff(v,2)/2+diff(u,2).*v/2;  % Frechet derivative
   Du=diff(u);                    % Needed to compute the BCs
   A.lbc=@(v) [v,diff(v)+Du(0)];
   A.rbc=@(v) diff(v)+Du(infty)-1;
   v = A\(-N(u));                 % Solve the linearized BVP      
   nrmv = norm(v)
   u = u + v;
end
du=diff(u);
d2u=diff(u,2);

fprintf('Kowarth reports f''''(0) = 0.332057.\n'); 
fprintf('Value computed here is f''''(0) = %7.10f.\n',d2u.vals(1)); 

要点：

1. Chebfun系统大量使用封装和重载技术；

2. 第三个层面的A算子就是第二个层面的J算子，区别在于：在第三个层面，A算子显式地写出，而在第二个层面上J通过对N算子使用diff，自动计算出。

3. 使用Chebfun求解Blasius方程，解对截断区间不敏感！！(比BVP4C好)

本文的计算基于chebfun_v4.2.2194，进一步研究参见使用Chebfun求解Blasius方程(二)。

参考文献：

1). Driscoll, Tobin A.; Bornemann, Folkmar; Trefethen, Lloyd N. (Nov. 22, 2008). "The Chebop system for automatic solution of differentialequations". BIT Num. Math. 48: 701–723.

2). A. Birkisson and T. A. Driscoll, AutomaticFrechet differentiation for the numerical solution of boundary-value problems,ACM Transactions on Mathematical Software.

3). RB Platte, LN Trefethen, Chebfun: A newkind of numerical computing.

4). http://www2.maths.ox.ac.uk/chebfun/guide/html/guide10.html

5). http://www2.maths.ox.ac.uk/chebfun/guide/html/guide7.shtml