《有限差分算例》由会员分享,可在线阅读,更多相关《有限差分算例(3页珍藏版)》请在金锄头文库上搜索。
1、有方程:c = 1, Re = 10 ;dudu1 d 2u+ c dtdx Re dx 2初始条件为:u (x,0)-罟(-1 x 0)(0 x 1)边界条件为:u (-1,t) 1 ,u(1,t) 0 ,给出定常解。一、算法描叙1. 网格剖分取x G -1,1, t G 0,5 (t值可以取大一些,到后期数值变化不大,可以取一定精度)x 方向 dx 0.1 , t 方向上 dt 0.001 。2. 差分格式iu j +1 - u jiidtdu )j uj ujI i+1i1 , (dx 丿2dxid 2u )j uj uj + ujii1,、dx 2 丿dx 2idtr1 2dxr2dt
2、10dx 2uj+1(r2r1)uj +(12*r2)uj +(r2+r1)uj 。 ii +1ii 13. 初边值处理直接给出。4. 稳定值判断当上下两层值之差的最大值不超过105,即认为结果达到要求。、程序:%*清理内存数据clear all;clc;/r/l/l *t7 a r I / V X 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1%*网格剖分*x0=-l;%空间起点xn=l;%空间终点t0=0;%时间起点tn=5;%时间终点(可以取得大一些) dx=0.1;%
3、空间步长 dt=0.001;%时间步长N=(xn-x0)/dx; T=(tn-t0)/dt;u=zeros(N+1,T+1);%设置整个网格上的0矩阵Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Txfor i=1:N/2+1u(i,1)=1;endfor i=N/2+2:N+1 u(i,1)=0;end %以上得到第1 层的值*1* *1* *1* *1* *1* *1* *1
4、* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Txfor j=1:T+1u
5、(1,j)=1; u(N+1,j)=0;end*1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx
6、Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Txx=x0:dx:xn; pc=plot(x,u(:,1),r,linewidth,2); axis(x0,xn,0,1);% 设置坐标区域 set(pc,EraseMode,xor);% 设置动画方案 xlabel(x);ylabel(u);title(u随时间变化曲线);添加标签Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx T
7、x Tx Tx Tx Tx Tx Txr1=dt/(dx*2); r2=dt/(10*dxA2); for j=1:T for i=2:N u(i,j+1)=(r2-r1)*u(i+1,j)+(1-2*r2)*u(i,j)+(r2+r1)*u(i-1,j);%由前层得到第2层及其以上各层数值endTx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Txe=abs(u(2:N,j+1)-u(2:N,j);e1=max(e);if e11e-5 & e1=0breakend*1* *1* *1*
8、 *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Tx Txset(pc,XData,x,YData,u(:,
9、j); drawnow;%刷新屏幕 pause(0.001)%刷新时间间隔end*1* *1*rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* rT* r
10、T* rT* rT* rT* rT* rT*fprintf(le-5精度数值解的稳定时刻t(s):) disp(j-1)*dt);U=u(l:N+l,j);fprintf(1e-5 精度数值解:) disp(U);结果1已-5精度數值解的稳定时刻土仪):2.42701-5 精度数值解: Columns 1 through 111.00001.00001.00001. 00001.00001.00001.00001.00001.00000. 99990.9999CoIujiltls 12 through 210.99980. 99950.99900. 99770.99450. 98560.96020. 88560.6637