函数文件:
1 function [x,n,flag]=sor(A,b,eps,M,max1) 2 %sor函数为用松弛迭代法求解线性方程组 3 %A为线性方程组的系数矩阵 4 %b为线性方程组的常数向量 5 %eps为精度要求 6 %M为超弛因子 7 %max1为最大迭代次数 8 %u为线性方程组的解 9 %n为迭代次数 10 %flag为指标变量,flag='OK!'表示迭代收敛达到指标要求 11 %flag='fail!'表示迭代失败 12 if nargin<5 13 max1=10000; 14 end 15 if nargin<4 16 M=1; 17 end 18 if nargin<3 19 eps=1e-11; 20 end 21 k=length(A); 22 n=0; 23 x=zeros(k,1); 24 y=zeros(k,1); 25 flag='OK!'; 26 while 1 27 y=x; 28 for i=1:k 29 z=b(i); 30 for j=1:k 31 if j~=i 32 z=z-A(i,j)*x(j); 33 end 34 end 35 if abs(A(i,i))<1e-10 | n==max1 36 flag='fail!'; 37 return; 38 end 39 z=z/A(i,i); 40 x(i)=(1-M)*x(i)+M*z; 41 end 42 if norm(y-x,inf)<eps 43 break; 44 end 45 n=n+1; 46 end
脚本文件:
1 tic; 2 clear 3 clc 4 N=100; 5 h=1/N; 6 x=0:h:1; 7 for i=1:length(x) 8 Accurate(i)=sin(4*pi*x(i)); 9 end 10 Accurate=Accurate';%精确解 11 A=diag(2/h^2*ones(N+1,1))+diag(-1/h^2*ones(N,1),1)+diag(-1/h^2*ones(N,1),-1); 12 A(1,1)=1; 13 A(1,2)=0; 14 A(N+1,N+1)=1; 15 A(N+1,N)=0; 16 b=zeros(N+1,1); 17 for i=1:N-1 18 b(i+1,1)=16*pi^2*sin(4*pi*x(i+1)); %右端函数 19 end 20 u0=zeros(N+1,1); 21 [u,n]=GaussSeid(A,b,u0) 22 numerial=u;%数值解 23 24 toc; 25 figure(1) 26 plot(x,Accurate,'r *',x,numerial,'g v'); 27 legend('Accurate','numerial'); 28 xlabel('x'); 29 ylabel('y'); 30 grid on; 31 toc; 32 figure(2) 33 plot(x,numerial-Accurate,'r *'); 34 legend('error'); 35 xlabel('x'); 36 ylabel('y'); 37 grid on;