• 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;
  • 相关阅读:
    [湖南集训]谈笑风生
    【SCOI2010】序列操作
    ●BZOJ 3994 [SDOI2015]约数个数和
    ●BZOJ 3309 DZY Loves Math
    ●UOJ 21 缩进优化
    ●BZOJ 2693 jzptab
    ●BZOJ 2154 Crash的数字表格
    ●BZOJ 3529 [Sdoi2014]数表
    ●2301 [HAOI2011] Problem b
    ●BZOJ 2820 YY的GCD
  • 原文地址:https://www.cnblogs.com/xtu-hudongdong/p/6507014.html
Copyright © 2020-2023  润新知