zoukankan      html  css  js  c++  java
  • 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;

  • 相关阅读:
    005. Asp.Net Routing与MVC 之三: 路由在MVC的使用
    004. Asp.Net Routing与MVC 之二: 请求如何激活Controller和Action
    001. Asp.Net Routing与MVC 之(基础知识):URL
    002. Asp.Net Routing与MVC 之(基础知识):HttpModule 与 HttpHandler
    003. Asp.Net Routing与MVC 之一: 请求如何到达MVC
    Factory
    decorator
    Java 单例真的写对了么?
    Dubbo Jackson序列化使用说明
    使用JavaConfig方式配置dubbox
  • 原文地址:https://www.cnblogs.com/xtu-hudongdong/p/6569040.html
Copyright © 2011-2022 走看看