使用Chebfun求解Blasius方程(二)

在上一篇文章使用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.