• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    迪恩网络公众号

Matlab:椭圆方程的导数边值问题

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

 

 

 

 

 

 

 

 

 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;

 


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
Delphi中匿名方法动态绑定事件发布时间:2022-07-18
下一篇:
delphi中Webbrowser疑难问题集锦转发布时间:2022-07-18
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap