抛物型差分（一维）

抛物型差分（一维）

从上到下依次为显式、隐式、精确解、六点差分得到的解

//Explicit.m

M = 10;

h = 1/( M + 1);

tao = 0.5 * h * h;

r = tao /( h * h);

I = [ 0 1 ];

T = [ 0 4 ];

init = zeros(1 , M + 2);

init(1) = sin(0);

init( M + 2) = sin( pi);

x = 2: M + 1;

x = x * h;

init( 2 : M + 1) = sin( pi * x) + x.*( 1 - x);

SOL = zeros(1 ,M + 2);

SOLPLOT = zeros( (T(2) - T(1))/tao ,M + 2);

for i = 1 : (T(2) - T(1))/tao

    SOL(1) = 0;

    SOL(M + 2) = 0;

    for j = 2 : M + 1

        SOL(j) = init(j - 1) * r + init(j) * (1 - 2*r) + init(j + 1) * r ...

            +  2 * tao;

    end

    init = SOL;

    SOLPLOT(i ,:) = init;

end

SOLPLOT = SOLPLOT';

time =  ones( M + 2 ,1) * linspace( T(1), T(2), (T(2) - T(1))/tao);

space = linspace( I(1), I(2) , M + 2)' * ones( 1 , (T(2) - T(1))/tao);

surf( time ,space ,SOLPLOT,'EdgeColor','none');

 

x1 = linspace( 0 ,1, M + 2)';

t1 = linspace( 0, 4 , (T(2) - T(1))/tao);

z = sin( pi * x1)*exp( -pi^2 * t1) + x1.*(1 - x1) * ones(1 ,length(t1));

figure;

mesh( t1, x1, z);

 

//Imp.m ，其余同Explicit.m

A = diag( ones(1, M) * ( 1 + 2 * r) ) + diag( ones( 1, M -1 ) * (-r),1)...

    + diag( ones(1, M -1) * (-r) , -1);

b = zeros( 1 ,M);

for i = 1 : (T(2) - T(1))/tao

    b = init( 2: M + 1);

    b(1) = b(1) + 0;

    b(M) = b(M) + 0;

    b = b + 2 * tao;

    init(2: M + 1) = A\b';

    init(1) = 0;

    init(M + 2) = 0;

    SOLPLOT(i, :) = init;

end

 

//Crank_Nicolson.m，其余同Explicit.m

theta = 1/2;

A = diag( ones(1, M) * ( 1 + 2 * r * theta) ) + diag( ones( 1, M -1 ) * (-r * theta),1)...

    + diag( ones(1, M -1) * (-r * theta) , -1);

b = zeros( 1 ,M);

for i = 1 : (T(2) - T(1))/tao

    SOL = init;

    for j = 2 : M + 1

        b( j - 1) =  theta * r * SOL(j-1) + (1 - 2 * theta * r) * SOL(j) +...

            theta * r * SOL(j + 1) + 2 * tao;

    end

    b(1) = b(1) + theta * r * 0;

    b(M) = b(M) + theta * r * 0;

    init(2: M + 1) = A\b';

    init(1) = 0;

    init(M + 2) = 0;

    SOLPLOT(i, :) = init;

end

 

注：1.现只做了一维的

   2.考虑的边界也是很特殊的，为长方的

   3.上述差分格式解稳定的条件是 r<= 1/2