- 1、本文档共7页,可阅读全部内容。
- 2、原创力文档(book118)网站文档一经付费(服务费),不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
- 3、本站所有内容均由合作方或网友上传,本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺!文档内容仅供研究参考,付费前请自行鉴别。如您付费,意味着您自己接受本站规则且自行承担风险,本站不退款、不进行额外附加服务;查看《如何避免下载的几个坑》。如果您已付费下载过本站文档,您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“版权申诉”(推荐),也可以打举报电话:400-050-0827(电话支持时间:9:00-18:30)。
查看更多
一维扩散方程有限差分法matlab
一维扩散方程的有限差分法——计算物理实验作业七陈万 物理学2013级 目:编程求解一维扩散方程的解取。输出t=1,2,...,10时刻的x和u(x),并与解析解u=exp(x+0.1t)作比较。主程序:% 一维扩散方程的有限差分法clear,clc;%定义初始常量a1 = 1; b1 = 1; c1 = 0; a2 = 1;b2 = -1; c2 = 0;a0 = 1.0; t_max = 10; D = 0.1; h = 0.1; tao = 0.1;%调用扩散方程子函数求解u = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2);子程序1:function output = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2)% 一维扩散方程的有限差分法,采用隐式六点差分格式(Crank-Nicolson)% a0: x的最大值% t:_max: t的最大值% h: 空间步长% tao: 时间步长% D:扩散系数% a1,b1,c1是(x=0)边界条件的系数;a2,b2,c2是(x=a0)边界条件的系数x = 0:h:a0;n = length(x);t = 0:tao:t_max;k = length(t); P = tao * D/h^2;P1 = 1/P + 1;P2 = 1/P - 1;u = zeros(k,n);%初始条件u(1,:) = exp(x);%求A矩阵的对角元素dd = zeros(1,n);d(1,1) = b1*P1+h*a1;d(2:(n-1),1) = 2*P1;d(n,1) = b2*P1+h*a2;%求A矩阵的对角元素下面一行元素ee = -ones(1,n-1);e(1,n-1) = -b2;%求A矩阵的对角元素上面一行元素ff = -ones(1,n-1);f(1,1) = -b1;R = zeros(k,n);%求R%追赶法求解for i = 2:k R(i,1) = (b1*P2-h*a1)*u(i-1,1)+b1*u(i-1,2)+2*h*c1;for j = 2:n-1 R(i,j) = u(i-1,j-1)+2*P2*u(i-1,j)+u(i-1,j+1);end R(i,n) = b2*u(i-1,n-1)+( b2*P2-h*a2)*u(i-1,n)+2*h*c2; M = chase(e,d,f,R(i,:)); u(i,:) = M'; plot(x,u(i,:)); axis([0 a0 0 t_max]); pause(0.1)endoutput = u;% 绘图比较解析解和有限差分解[X,T] = meshgrid(x,t);Z = exp(X+0.1*T);surf(X,T,Z),xlabel('x'),ylabel('t'),zlabel('u'),title('解析解');figuresurf(X,T,u),xlabel('x'),ylabel('t'),zlabel('u'),title('有限差分解');子程序2:function M = chase(a,b,c,f)% 追赶法求解三对角矩阵方程,Ax=f% a是对角线下边一行的元素% b是对角线元素% c是对角线上边一行的元素% M是求得的结果,以列向量形式保存n = length(b);beta = ones(1,n-1);y = ones(1,n);M = ones(n,1);for i = (n-1):(-1):1 a(i+1) = a(i);end% 将a矩阵和n对应beta(1) = c(1)/b(1);for i = 2:(n-1) beta(i) = c(i)/( b(i)-a(i)*beta(i-1) );endy(1) = f(1)/b(1);for i = 2:n y(i) = (f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1));endM(n) = y(n);for i = (n-1):(-1):1 M(i) = y(i)-beta(i)*M(i+1);endend结果:对比分析两图,结果令人满意。取出t_max时刻的u值分析:()有限差分解如下:2.7182818283.0041660243.3201169233.6692966684.0551999674.4816890704.9530324245.4739473926.0496474646.6858944427.389056099解析解如下:2.7145201222.99951097
您可能关注的文档
- 【三年高考两年模拟】2017年高考物理新课标一轮复习第五章 万有引力与航天 1讲 万有引力定律及天体运动.pptx
- 【下雨啦】201带电粒子在磁场中运动轨迹.ppt
- 【中考冲刺】2015年中考物理冲刺复习课(重难点突破+剖析重点实验):第五章 透视及其应用.ppt
- 【中考备战策略】2017年中考物理(人教版)总复习第19讲 电率(一).ppt
- 【中考突破】人教版2016年初中物理中考复习课件:《第章 光现象》ppt课件.ppt
- 【中考备战策】2017年中考物理(人教版)总复习第5讲 光现象.ppt
- 【中考备战策略】2017年中考物理(人教版)总复习阶段练(一).ppt
- 【人教版】高中化学选修三:1.1.1《能层与能、构造原理》ppt课件.ppt
- 【人教版】高中化学选修三:1.1.2《基态原子的外电子排布、原子轨道》ppt课件.ppt
- 【人教版】高化学选修3知识点总结.docx
文档评论(0)