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

Matlab:CrankNicolson方法求解线性抛物方程

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

 1 tic;
 2 clear
 3 clc
 4 M=[10,20,40,80,160,320,640];%x的步数
 5 K=M; %时间t的步数
 6 for p=1:length(M)
 7 hx=1/M(p);
 8 ht=1/K(p);
 9 r=ht/hx^2; %网格比
10 x=0:hx:1;
11 t=0:ht:1;
12 numerical=zeros(M(p)+1,K(p)+1);
13 numerical(:,1)=exp(x); %初始值
14 numerical(1,:)=exp(t); %边值
15 numerical(M(p)+1,:)=exp(t+1); %边值
16 a=-r/2*ones(M(p)-2,1);b=(1+r)*ones(M(p)-1,1);c=-r/2*ones(M(p)-2,1);
17 fun1=inline('exp(x+t)','x','t');
18 for i=1:length(x)
19     for j=1:length(t)
20 Accurate(i,j)=fun1(x(i),t(j));
21     end
22 end
23 d=r/2*ones(M(p)-2,1);e=(1-r)*ones(M(p)-1,1);f=r/2*ones(M(p)-2,1);
24 B=diag(d,-1)+diag(e,0)+diag(f,1);
25 fun2=inline('0','x','t');
26 for i=1:M(p)-1
27 for k=1:K(p)
28     f(i,k)=ht*fun2(x(i+1),t(k)+ht/2);
29 end
30 end
31 for k=1:K(p)
32 f(1,k)=r/2*(numerical(1,k+1)+numerical(1,k));
33 f(M(p)-1,k)=r/2*(numerical(M(p)+1,k+1)+numerical(M(p)+1,k));
34 end
35 for k=1:K(p)
36     right_vector=f(:,k)+B*numerical(2:M(p),k);
37     numerical(2:M(p),k+1)=chase(a,b,c,right_vector);
38 end
39 error=numerical(2:M(p),2:K(p))'-Accurate(2:M(p),2:K(p))';
40 error_inf(p)=max(max(error));
41 figure(p)
42 [X,Y]=meshgrid(x,t);
43 subplot(1,3,1)
44 mesh(X,Y,Accurate');
45 xlabel('x'),ylabel('t');zlabel('Accurate');
46 title('the image of Accurate result');
47 grid on
48 subplot(1,3,2)
49 mesh(X,Y,numerical');
50 xlabel('x'),ylabel('t');zlabel('numerical');
51 title('the image of numerical result');
52 grid on
53 subplot(1,3,3)
54 mesh(X,Y,numerical'-Accurate');
55 xlabel('x'),ylabel('t');zlabel('error');
56 title('the image of error result');
57 grid on
58 end
59 for k=2:length(M)
60     H=error_inf(p-1)/error_inf(p);
61 E_inf(k-1)=log2(H);
62 end
63 figure(length(M)+1)
64 plot(1:length(M)-1,E_inf,'-r v');
65 ylabel('误差阶数');
66 title('Crank nicolson 误差阶数');
67 grid on
68 toc;

 


鲜花

握手

雷人

路过

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

请发表评论

全部评论

专题导读
上一篇:
Octave/Matlab初步学习发布时间:2022-07-22
下一篇:
Matlabtxt内容替换函数fgetlfseek发布时间:2022-07-22
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

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

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

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