一维抛物线偏微分方程数值解法(4)(附图及matlab程序)

1、一维抛物线偏微分方程数值解法（4）上一篇参看 一维抛物线偏微分方程数值解法（3）（附图及matlab程序）解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）Ut-Uxx=0, 0 （ Ut-aUxx=f(x,t,a>0 U(x,0=ex, 0<=x<=1,U(0,t=et,U(1,t=e(1+t, 0 精确解为：U(x,t=e(x+t;用紧差分格式：此种方法精度为o(h12+h24,无条件差分稳定；一：用追赶法解线性方程组（还可以用迭代法解）Matlab程序为：function u p e x t=JCHGS(h1,h2,m,n%紧差分格式解一维抛物线型偏微分方

2、程%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m,n分别为空间，时间网格数%p为精确解，u为数值解,e为误差x=(0:m*h1+0; x0=(0:m*h1;%定义x0,t0是为了f(x,t=0的情况%t=(0:n*h2+0; t0=(0:n*h2+1/2*h2;syms f;for(i=1:n+1for(j=1:m+1f(i,j=0; %f(i,j=f(x0(j,t0(i=0%endendfor(i=1:n+1u(i,1=exp(t(i;u(i,m+1=exp(1+t(i;endfor(i=1:m+1u(1,i=exp(x(i;endr=h2/(h1*h1;for(i=1

3、:n %外循环，先固定每一时间层，每一时间层上解一线性方程组%a(1=0;b(1=5/6+r;c(1=1/12-r/2;d(1=(r/2-1/12*u(i+1,1+.(1/12+r/2*u(i,1+(5/6-r*u(i,2+(1/12+r/2*u(i,3+.h2/12*(f(i,1+10*f(i,2+f(i,3;for(k=2:m-2a(k=1/12-r/2;b(k=5/6+r;c(k=1/12-r/2;d(k=h2/12*(f(i,k+.10*f(i,k+1+f(i,k+2+(1/12+r/2*(u(i,k+u(i,k+2+(5/6-r.*u(i,k+1;%输入部分系数矩阵，为0的矩阵元素不

4、输入%一定要注意输入元素的正确性enda(m-1=1/12-r/2;b(m-1=5/6+r;d(m-1=(1/12+r/2*(u(i,m-1+u(i,m+1+.(5/6-r*u(i,m+(r/2-1/12*u(i+1,m+1+ .h2/12*(f(i,m-1+10*f(i,m+f(i,m+1;for(k=1:m-2 %开始解线性方程组 消元过程a(k+1=-a(k+1/b(k;b(k+1=b(k+1+a(k+1*c(k;d(k+1=d(k+1+a(k+1*d(k;endu(i+1,m=d(m-1/b(m-1; %回代过程% for(k=m-2:-1:1u(i+1,k+1=(d(k-c(k*u(

5、i+1,k+2/b(k;endendfor(i=1:n+1for(j=1:m+1p(i,j=exp(x(j+t(i; %p为精确解e(i,j=abs(u(i,j-p(i,j;%e为误差endendu p e x t=JCHGS(0.1,0.005,10,200; surf(x,t,e>> title('误差'运行约43秒；u p e x t=JCHGS(0.1,0.01,10,100;surf(x,t,e 20多秒；u p e x t=JCHGS(0.2,0.04,5,25;surf(x,t,e 3秒；此方法精度很高；二：g-s迭代法求解线性方程组Matlab程序f

6、unction u e p x t k=JCFGS1(h1,h2,m,n,kmax,ep% 解抛物线型一维方程 格式 （Ut-aUxx=f(x,t,a>0%用g-s(高斯-赛德尔迭代法解%kmax为最大迭代次数%m,n为x,t方向的网格数，例如（2-0）/0.01=200;%e为误差，p为精确解syms temp;u=zeros(n+1,m+1;x=0+(0:m*h1;t=0+(0:n*h2;for(i=1:n+1u(i,1=exp(t(i;u(i,m+1=exp(1+t(i;end for(i=1:m+1u(1,i=exp(x(i;endfor(i=1:n+1for(j=1:m+1f(

7、i,j=0;endenda=zeros(n,m-1;r=h2/(h1*h1; %此处r=a*h2/(h1*h1）；a=1for(k=1:kmaxfor(i=1:nfor(j=2:mtemp=(1/12+r/2*(u(i,j-1+u(i,j+1+(5/6-r*u(i,j+.h2/12*(f(i,j-1+10*f(i,j+f(i,j+1+(r/2-1/12*(u(i+1,.j-1+u(i+1,j+1/(5/6+r;a(i+1,j=(temp-u(i+1,j*(temp-u(i+1,j;u(i+1,j=temp;%此处注意是u(i+1,j,而不是u(i+1,j+1%endenda(i+1,j=sqrt(a(i+1,j;if(k>kmaxbreak;endif(max(max(a break; endendfor(i=1:n+1for(j=1:m+1p(i,j=exp(x(j+t(i;e(i,j=abs(u(i,j-p(i,j;endendu e p x t k=JCFGS1(0.1,0.