Demo entry 6737220

MachineExperiment6_jacobi_solution

   

Submitted by Jone on Apr 30, 2018 at 04:57
Language: C++. Code size: 1.9 kB.

function [y,k]=(eps,h,delta)  %a=0.5 (a,eps,h,delta)
% =============  Initialization:
n=1.0/h;     %步长h
A=ones(n-1);   %矩阵A
y=zeros(n-1,1);   %JacobiY
z=zeros(n-1,1);   %临时z
k=0;           %迭代次数

% =============  聪明型画矩阵A:
A=diag(-(2*eps+h)*ones(1,n))+diag(eps*ones(1,n-1),-1)+diag((eps+h)*ones(1,n-1),1);
% =============  傻瓜型画矩阵A:
% for i=1:n-1
%     for j=1:n-1
%         A(i,j)=0;
%     end
% end
% 
% for i=1:n-1
%     A(i,i)=-(2*eps+h);
% end
% 
% for i=1:n-1
%     for j=1:n-1
%         if i==j+1
%             A(i,j)=eps;
%         end
%         if i==j-1
%             A(i,j)=eps+h;
%         end
%     end
% end
% =============  系数矩阵b :
b=zeros(n-1,1);
for i=1:n-2
    b(i,1)=0.5*h^2;
end
b(n-1,1)=0.5*h^2-eps-h;

% =============  构造对角矩阵D :
D=zeros(n-1);
for i=1:n-1
    D(i,i)=A(i,i);
end

% =============  计算A的下三角矩阵L :
L=zeros(n-1);
for i=1:n-1
    for j=1:n-1
if i>j
    L(i,j)=-A(i,j);
end
    end
end

% =============  计算A的上三角矩阵U :
U=zeros(n-1);
for i=1:n-1
    for j=1:n-1
        if i<j
            U(i,j)=-A(i,j);
        end
    end
end

% =============  Jacobi核心程序 :
B=D\(L+U);           %Jacobi迭代法的迭代矩阵 B=D^-1(L+U)
g=D\b;               %g=D^-1*b
while 1
    z=B*y+g;         %临时z
    if norm(z-y,inf)<delta   % norm(z-y,inf),相当于返回max(abs(z-y)) 
           break;            % norm(A,'inf'),返回矩阵的无穷泛数,也就是最大一行的和
    end
    y=z;           
    k=k+1;
end
k                  %输出显示迭代次数k

% ============= Drawing an image :
x=linspace(0,1);      %[0,1]区间默认100等分
truey=(1-0.5)/(1-exp(-1/eps))*(1-exp(-x./eps))+x.*0.5;   %真实解,a=0.5
%e=truey+y;
figure;               %打开绘画窗口
plot(100*x,truey,'g','LineWidth',2.0);      
hold on;               %保持曲线
plot(y,'r','LineWidth',1.0);        %绘制解图像,颜色绿,线宽1.0
legend('True Y','Jacobi Y');              %图例
grid;                 %开网格

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).