导数边界

Posted 2020-09-05 胡冬冬

 

 

 

 

 

 

 

 

 1 tic;
 2 clear
 3 clc
 4 N=8;
 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/h1*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/h1*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/h1*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/h1*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))+1/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     end
55 end
56 error=Numerical-numerical;
57    numerical=Numerical;
58 end
59 for i=1:length(x)
60     for j=1:length(y)
61    Accurate(i,j)=fun(x(i),y(j));
62     end
63 end
64 Error=Accurate\'-Numerical\';
65 [X,Y]=meshgrid(x,y);
66 subplot(1,3,1)
67 mesh(X,Y,Accurate\');
68 xlabel(\'x\');ylabel(\'y\');zlabel(\'Accurate\');
69 grid on
70 subplot(1,3,2)
71 mesh(X,Y,Numerical\');
72 xlabel(\'x\');ylabel(\'y\');zlabel(\'Numerical\');
73 grid on
74 subplot(1,3,3)
75 mesh(X,Y,Error);
76 xlabel(\'x\');ylabel(\'y\');zlabel(\'error\');
77 grid on
78 toc;