资源描述:
《偏微分方程数值解上机.doc》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、个人收集整理勿做商业用途偏微分方程数值解上机实习数值求解二维扩散方程的初边值问题古典显式格式:将原格式化为:附源程序:%——---—-—--——---—--—-—-———-———-----—---——-—-运用古典显式差分格式求解二维扩散方程的初边值问题;functiongdxs(ti,h,t)%--—-——-——————------———---—-—-------—-——-———ti:时间步长;%-——-——-—-—-————---———-——-———------——-—---——h:空间步长;k=t/ti;m=1/h+1;r=ti/h^2;%———————-—---——-
2、——-—-—--—----——r为网格比;w=ones(m,m);u=ones(m,m);fori=2:m-1forj=2:m—1u(i,j)=sin(pi*(i—1)*h)*sin(2*pi*h*(j-1));endendticforl=1:kfori=2:m—1forj=2:m-1w(i,j)=r*u(i-1,j)+r*u(i,j-1)+r*u(i+1,j)+r*u(i,j+1)+(1—4*r)*u(i,j);endendu=w;endtoct=tocumesh(u)个人收集整理勿做商业用途交替方向隐式格式(P-R格式):将原差分格式化为:代入边界条件,转化为三对角矩阵个
3、人收集整理勿做商业用途附追赶法源程序:%-----—-—————-—--—-—-—---—-——-—-—-————---—-—追赶法求解三对角方程组;functionx=zg(a,b,c,d)%--—-———-———-—-——-—--—————-—-—-—---——--——————a:方程组系数矩阵A的下对角元素;%—-—-—-——---——-—-——--—--—-——-—-———-——-—--—-——b:方程组系数矩阵A的主对角元素;%-—---———-———-————--——-—-—--——-——-—-——-----——c:方程组系数矩阵A的上对角元素;%----————
4、--—--——-——————-—--——--—---——----——-—d:追赶法所求方程的右端向量;%———--——---—--————-—-——--——--——-—--———----——-l:系数矩阵A所分解成的下三角阵L中的下对角元素了l(i);%—----—----—--—--—-————-—-—--—-——-----—-—-——-u:系数矩阵A所分解成的下三角阵U中的主对角元素了u(i);n=length(b);u(1)=b(1);y(1)=d(1);fori=1:n-1%-—-—--———-—------—-—-——---追赶法求解之追过程求解Ly=d;l(i)=
5、a(i)/u(i);u(i+1)=b(i+1)—l(i)*c(i);y(i+1)=d(i+1)-l(i)*y(i);endx(n)=y(n)/u(n);%————---———---—-—---——---追赶法求解之赶过程求解Uz=y;forj=n—1:—1:1ifu(j)==0break;elsex(j)=(y(j)—c(j)*x(j+1))/u(j);个人收集整理勿做商业用途endend%--—-—--————-———--—-——-—-—-----——-—-------—-——--运用P—R差分格式求解二维扩散方程的初边值问题;functionpr(ti,h,t)%-——-
6、--———--——--—-—-———-——---—-—--——--———---ti:时间步长h:空间步长;k=t/ti+1;m=1/h+1;r=ti/h^2;%-————-—--—--——--——-----—--——--r为网格比;w=ones(m,m);u=ones(m,m);%--——--—————--———-————-——输入初始值v=ones(m,m);fori=2:m-1forj=2:m-1u(i,j)=sin(pi*(i—1)*h)*sin(2*pi*h*(j-1));endend%-—--———----————---—--———输入用P—R差分格式求解的三对角矩
7、阵b=ones(1,m—2)*(2+2*r);a=-r*ones(1,m-3);c=-r*ones(1,m-3);A=zeros(m-2,m—2);fori=1:m—2A(i,i)=2-2*r;endfori=1:m-3A(i,i+1)=r;A(i+1,i)=r;endp=zeros(m-2,1);p(1)=2*r;p(m—2)=2*r;ticforl=1:kfori=2:m-1d1=A*u(i,2:m—1)’+p;d1=d1’;w(2:m—1,i)=zg(a,b,c,d1);%--————-——-——