导数边界

 1 tic;
 2 clear
 3 clc
 4 N=4;
 5 M=2*N;
 6 h1=2/M;
 7 h2=1/N;
 8 x=0:h1:2;
 9 y=0:h2:1;
10 fun=inline(‘exp(x)*sin(pi*y)‘,‘x‘,‘y‘);
11 f=inline(‘(pi^2-1)*exp(x)*sin(pi*y)‘,‘x‘,‘y‘);
12 lamda1=inline(‘1‘,‘y‘);
13 lamda2=inline(‘2*y‘,‘y‘);
14 lamda3=inline(‘2*x‘,‘x‘);
15 lamda4=inline(‘x^2‘,‘x‘);
16 kesai1=inline(‘0‘,‘y‘);
17 kesai2=inline(‘exp(2)*(1+2*y)*sin(pi*y)‘,‘y‘);
18 kesai3=inline(‘-pi*exp(x)‘,‘x‘);
19 kesai4=inline(‘-pi*exp(x)‘,‘x‘);
20 numerical=zeros(M+1,N+1);
21 Numerical=numerical;
22 error=eye(M+1,N+1);
23 while norm(error,inf) >= 1e-5
24 for j=1:N+1
25     for i=1:M+1
26         if i==1 & j==1
27             Numerical(i,j)=(f(x(i),y(j))+2/h1^2*numerical(i+1,j)+2/h2^2*numerical(i,j+1)+2/h1*kesai1(y(j))+2/h2*kesai3(x(i)))...
28                 /(2/h1^2+2/h2^2+2/h2*lamda1(y(j))+2/h2*lamda3(x(i)));%U(0,0)
29         elseif i==M+1 & j==1
30              Numerical(i,j)=(f(x(i),y(j))+2/h1^2*numerical(i-1,j)+2/h2^2*numerical(i,j+1)+2/h1*kesai2(y(j))+2/h2*kesai3(x(i)))...
31                 /(2/h1^2+2/h2^2+2/h2*lamda2(y(j))+2/h2*lamda3(x(i)));%U(m,0)
32         elseif i==1 & j==N+1
33             Numerical(i,j)=(f(x(i),y(j))+2/h1^2*numerical(i+1,j)+2/h2^2*numerical(i,j-1)+2/h1*kesai1(y(j))+2/h2*kesai4(x(i)))...
34                 /(2/h1^2+2/h2^2+2/h2*lamda1(y(j))+2/h2*lamda4(x(i)));%U(0,n)
35         elseif i==M+1 & j==N+1
36             Numerical(i,j)=(f(x(i),y(j))+2/h1^2*numerical(i-1,j)+2/h2^2*numerical(i,j-1)+2/h1*kesai2(y(j))+2/h2*kesai4(x(i)))...
37                 /(2/h1^2+2/h2^2+2/h2*lamda2(y(j))+2/h2*lamda4(x(i)));%U(m,n)
38         elseif i==1 & j>=2 & j<=N  %  0j
39             Numerical(i,j)=(f(x(i),y(j))+2/h1^2*numerical(i+1,j)+1/h2^2*numerical(i,j-1)+1/h2^2*numerical(i,j+1)+2/h1*kesai1(y(j)))...
40                 /(2/h1^2+2/h2^2+2/h1*lamda1(y(j)));
41         elseif j==1 & i>=2 & i<=M  % i0
42              Numerical(i,j)=(f(x(i),y(j))+2/h1^2*numerical(i-1,j)+1/h1^2*numerical(i+1,j)+2/h2^2*numerical(i,j+1)+2/h2*kesai3(x(i)))...
43                 /(2/h1^2+2/h2^2+2/h2*lamda3(x(i)));
44         elseif i==M+1 & j>=2 & j<=N  %  mj
45              Numerical(i,j)=(f(x(i),y(j))+2/h1^2*numerical(i-1,j)+1/h2^2*numerical(i,j-1)+1/h2^2*numerical(i,j+1)+2/h1*kesai2(y(j)))...
46                 /(2/h1^2+2/h2^2+2/h1*lamda2(y(j)));
47         elseif j==N+1 & i>=2 & i<=M  %  in
48              Numerical(i,j)=(f(x(i),y(j))+1/h1^2*numerical(i-1,j)+1/h1^2*numerical(i+1,j)+2/h2^2*numerical(i,j-1)+2/h2*kesai4(x(i)))...
49                 /(2/h1^2+2/h2^2+2/h2*lamda4(x(i)));
50         else
51              Numerical(i,j)=(f(x(i),y(j))+1/h1^2*numerical(i-1,j)+1/h2^2*numerical(i,j-1)+1/h1^2*numerical(i+1,j)+1/h2^2*numerical(i,j+1))...
52                  /(2/h1^2+2/h2^2);
53         end
54
55     end
56 end
57 error=Numerical-numerical;
58    numerical=Numerical;
59 end
60 for i=1:length(x)
61     for j=1:length(y)
62    Accurate(i,j)=fun(x(i),y(j));
63     end
64 end
65 Error=Accurate‘-Numerical‘;
66 [X,Y]=meshgrid(x,y);
67 mesh(X,Y,Error);
68 toc;