compdiff.m

来自「various matlab programs to slove various」· M 代码 · 共 25 行

function [u, x] = compdiff(alpha,beta,gamma,a,b,n,f)
h=(b-a)/(n+1); x=[a:h:b]; fx = eval(f);
A   = eye(n+2)+alpha*diag(ones(n+1,1),1)+alpha*diag(ones(n+1,1),-1);
rhs = 0.5*beta/h*(fx(4:n+1)-fx(2:n-1))+0.25*gamma/h*(fx(5:n+2)-fx(1:n-2));
if gamma == 0
   rhs=[0.5*beta/h*(fx(3)-fx(1)), rhs, 0.5*beta/h*(fx(n+2)-fx(n))];
   A(1,1:n+2) = zeros(1,n+2); 
   A(1,1) = 1; A(1,2)=alpha; A(1,n+1)=alpha; 
   rhs=[0.5*beta/h*(fx(2)-fx(n+1)), rhs];
   A(n+2,1:n+2) = zeros(1,n+2);
   A(n+2,n+2) = 1; A(n+2,n+1)=alpha; A(n+2,2)=alpha;
   rhs=[rhs, 0.5*beta/h*(fx(2)-fx(n+1))];
else
   rhs=[0.5*beta/h*(fx(3)-fx(1))+0.25*gamma/h*(fx(4)-fx(n+1)), rhs];
   A(1,1:n+2) = zeros(1,n+2); 
   A(1,1) = 1; A(1,2)=alpha; A(1,n+1)=alpha;
   rhs=[0.5*beta/h*(fx(2)-fx(n+1))+0.25*gamma/h*(fx(3)-fx(n)), rhs];
   rhs=[rhs,0.5*beta/h*(fx(n+2)-fx(n))+0.25*gamma/h*(fx(2)-fx(n-1))];
   A(n+2,1:n+2) = zeros(1,n+2);
   A(n+2,n+2) = 1; A(n+2,n+1)=alpha; A(n+2,2)=alpha;
   rhs=[rhs,0.5*beta/h*(fx(2)-fx(n+1))+0.25*gamma/h*(fx(3)-fx(n))];
end
u = A \ rhs';
return