在上一篇文章使用Chebfun求解Blasius方程(一)里，我们使用chebfun求解了Blasius方程。 由于Blasius方程定义在半无界区间，因此我们将区间进行截断以求解，也就是说，在无穷远处的边界条件f’(+∞)=1被f’(infty)=1替代，此处infty是一个比较大的数字，如10，20，100等，但是并非是无穷大。
Chebfun能很好地表示具有可去奇点的函数，利用这一特点，我们试图对区间不截断而整体求解。基本思想是，将原方程的解变换到一个新的坐标系，而在这个新坐标系中求解区间从半无限区间变为有限区间, 并且方程的解在这个有限区间上没有奇性。当我们在这个“合适的”坐标系求解出变换后方程的解后，再变换到原坐标系而得到原方程的解。注意，Blasius方程的解有两个奇点：第一，方程的解定义在半无界区间上；第二，方程的解f(η)当η趋于正无穷大时是无穷大(f(η)→η)。
在变换
之下，原Blasius方程变为
若已知u(x),则原方程的解可以表示为
使用Chebfun求解关于u(x)的方程如下：
clc;
clear;
cheboppref('display','iter');
N=chebop(0,1);
x=chebfun('x',[0,1]);
N.op=@(u) -2*x.^4.*diff(u,3)+x.^2.*u.*diff(u,2)...% equation of u(x)
+(x-x.^2-12*x.^3).*diff(u,2)+2*x.*u.*diff(u)+2*(1-x-6*x.^2).*diff(u);
N.lbc=@(u) diff(u);
N.rbc=@(u) [u diff(u)-1];
N.init = -0.5+0.5*x.^2;
u=N\0;
figure(1)
plot(u);
title('Blasius equation in another coordinate ')
xlabel('x')
ylabel('u(x)')
infty=20;
eta=chebfun('eta',[0,infty]);
f=u(1./(eta+1))+eta; % coordinate transformation
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));
figure(2)
subplot(1,2,1);
plot(f);
title('Blasius equation')
xlabel('\eta')
ylabel('f(\eta)')
subplot(1,2,2);
plot(df);
axis([0 infty 0 1.1]);
title('Blasius equation')
xlabel('\eta')
ylabel('df/d\eta')
shg注意:
1. 迭代6步达到收敛精度；
2. 初始猜测解是多项式；
3.在chebfun系统上，从u(x)得到f(η)的转变，也就是复合函数的表示；
计算结果如下：
Iter. || update || length(update) stepsize length(u)
---------------------------------------------------------------------------
1 6.105e-001 129 1 129
2 1.620e-001 146 1 146
3 1.755e-002 164 1 164
4 1.980e-004 140 1 164
5 2.323e-008 100 1 164
6 1.552e-015 2 1 164
6 iterations
Final residual norm: 2.67e-006 (interior) and 8.74e-012 (boundary conditions).
Kowarth reports f''(0) = 0.332057.
Value computed here is f''(0) = 0.3320573362.
