Matlab:椭圆方程的导数边值问题
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;
胡冬冬
