Matlab-9:中心差分方法解常微分算例（SOR完整版）

函数文件：

 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;