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

Matlab:五点差分方法求解椭圆方程非导数边值问题

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

 

 

 

 

 差分格式脚本文件:

 1 tic;
 2 clear
 3 clc
 4 M=32;%x的步数
 5 N=16;%y的步数
 6 h1=1/M;%x的步长
 7 h2=1/N;%y的步长
 8 x=0:h1:1;
 9 y=0:h2:1;
10 u=zeros(M+1,N+1);%给数值解分配内存单元
11 U=u;%给精确解分配内存单元
12 u(1,:)=y.^3;%y边值
13 u(M+1,:)=1+y.^3;%y边值
14 u(:,1)=x.^3;%x边值。uo
15 u(:,N+1)=1+x.^3;%x边值。un
16 for i=1:M+1
17     for j=1:N+1
18     Accurate(i,j)=x(i)^3+y(j)^3;%精确解
19     end 
20 end
21 fun=inline(\'-6*(x+y)\',\'x\',\'y\');
22 for i=1:M-1
23 for j=1:N-1
24         f(i,j)=fun(x(i+1),y(j+1));
25 end
26 end
27 Numerical=u;
28 error=eye(M+1,N+1);
29  while norm(error,inf) >= 1e-5
30 for i=2:M
31     for j=2:N
32         Numerical(i,j)=(f(i-1,j-1)+N^2*u(i,j-1)+M^2*u(i-1,j)+M^2*u(i+1,j)+N^2*u(i,j+1))/(2*M^2+2*N^2);    
33     end
34 end
35  error=Numerical-u;
36    u=Numerical;
37  end
38 [X,Y]=meshgrid(x,y);
39 subplot(1,3,1)
40 mesh(X,Y,Numerical\');
41 title(\'the image of Numerical solution\')
42 xlabel(\'x\');ylabel(\'y\');zlabel(\'u\');
43 subplot(1,3,2)
44 mesh(X,Y,Accurate\');
45 title(\'the image of Accurate solution\')
46 xlabel(\'x\');ylabel(\'y\');zlabel(\'U\');
47 subplot(1,3,3)
48 mesh(X,Y,(Numerical-Accurate)\');
49 title(\'the image of error solution\')
50 xlabel(\'x\');ylabel(\'y\');zlabel(\'error\');
51 toc;

 

效果图:

紧差分格式:

 1 tic;
 2 clear
 3 clc
 4 M=100;%x的步数
 5 N=100;%y的步数
 6 h1=1/M;%x的步长
 7 h2=1/N;%y的步长
 8 x=0:h1:1;
 9 y=0:h2:1;
10 u=zeros(M+1,N+1);%给数值解分配内存单元
11 U=u;%给精确解分配内存单元
12 u(1,:)=y.^3;%y边值
13 u(M+1,:)=1+y.^3;%y边值
14 u(:,1)=x.^3;%x边值。uo
15 u(:,N+1)=1+x.^3;%x边值。un
16 for i=1:M+1
17     for j=1:N+1
18     Accurate(i,j)=x(i)^3+y(j)^3;%精确解
19     end 
20 end
21 b1=5/3*(M^2+N^2);
22 b2=1/6*(5*M^2-N^2);
23 b3=1/6*(5*N^2-M^2);
24 f=inline(\'-6*(x+y)\',\'x\',\'y\');
25 for i=1:M-1
26 for j=1:N-1
27         ABf(i,j)=(1/144)*(f(x(i),y(j))+10*f(x(i),y(j+1))+f(x(i),y(j+2))...
28             +10*f(x(i+1),y(j))+100*f(x(i+1),y(j+1))+10*f(x(i+1),y(j+2))...
29             +f(x(i+2),y(j))+10*f(x(i+2),y(j+1))+f(x(i+2),y(j+2)));
30 end
31 end
32 Numerical=u;
33 error=eye(M+1,N+1);
34  while norm(error,inf) >= 1e-10
35 for i=2:M
36     for j=2:N
37         Numerical(i,j)=(ABf(i-1,j-1)+1/20*b1*Numerical(i-1,j-1)+b3*Numerical(i,j-1)+1/20*b1*Numerical(i+1,j-1)...
38             +b2*Numerical(i-1,j)+b2*u(i+1,j)...
39             +1/20*b1*u(i-1,j+1)+b3*u(i,j+1)+1/20*b1*u(i+1,j+1))/b1;    
40     end
41 end
42  error=Numerical-u;
43    u=Numerical;
44  end
45 [X,Y]=meshgrid(x,y);
46 subplot(1,3,1)
47 mesh(X,Y,Numerical\');
48 title(\'the image of Numerical solution\')
49 xlabel(\'x\');ylabel(\'y\');zlabel(\'u\');
50 subplot(1,3,2)
51 mesh(X,Y,Accurate\');
52 title(\'the image of Accurate solution\')
53 xlabel(\'x\');ylabel(\'y\');zlabel(\'U\');
54 subplot(1,3,3)
55 mesh(X,Y,(Numerical-Accurate)\');
56 title(\'the image of error solution\')
57 xlabel(\'x\');ylabel(\'y\');zlabel(\'error\');
58 toc;

效果图:

 


鲜花

握手

雷人

路过

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

请发表评论

全部评论

专题导读
上一篇:
多个接口出现同名函数的处理(C#,Delphi)发布时间:2022-07-18
下一篇:
Delphi中关于Move的操作发布时间: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