使用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好)

补充：综合编程 , 其他综合 ,