3

私は matlab でクランク ニコルソン法を実装しようとしていますが、境界条件 (つまり、u(0,t)=u(N,t)=0) なしで実装を機能させることができました。私が抱えている問題は、境界条件を追加することです。私の現在の実装では、境界条件が考慮されていないようです。

ここに私の現在の実装があります: CN メソッド:

function [ x, t, U ] = Crank_Nicolson( vString, fString, a, N, M,g1,g2 )
%The Crank Nicolson provides a solution to the parabolic equation provided
%   The Crank Nicolson method uses linear system of equations to solve the
%   parabolic equation.

%Prepare the grid and grid spacing variables. 
dt = 1/M;
t = dt * [0:M];
h = 1/N;
x = 2 + h * [0:N]';%Shift x by 2 that way we have 2 <= x <= 3

%Prepare the matrix that will store the solutions over the grid
U = zeros(N+1, M+1);
%This will fill the first column with the initial condition. feval will
%evaluate the initial condition at all values of x
U(:,1) = feval(vString, x);
%This fills the boundary conditions. One boundary condition goes on the
%first row the other boundary condition goes on the last row
U(1,:) = feval(g1, t);
U(end,:) = feval(g2, t);

%The loop that will populate the matrix with the solution
n = 1:N+1;%Start at 2 since n=1 is the initial condition
e = ones(N+1,1);
B = spdiags([-1*e 2*e -1*e],-1:1, N+1, N+1)*(1/h^2);
A = (speye(N+1)+((a*dt)/2)*B);
X = (speye(N+1)-((a*dt)/2)*B);
R = chol(A);%Choleski decomposition
for m=2:M+1
    %The linear system is solved. 
    b = X*U(n,m-1) + dt * feval(fString, x(n), (t(m)+t(m-1))*0.5);
    b = R'\b;
    U(n,m) = R\b;
end
end

この実装は、境界条件が問題にならない場合に機能することを知っています。足りないものはありますか?また、私は matlab に比較的慣れていないため、一般的な matlab 形式の提案があれば喜んでお知らせします。

プロジェクト全体に興味がある場合は、これをダウンロードしてください

4

0 に答える 0