维抛物线方程数值解法(ADI隐式交替法)方法

PAGE/NUMPAGESADI隐式交替法三种解法及误差 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 （一般的教材上只说第一种）理论部分参看孙志忠：偏微分方程数值解法注意：最好不要直接看程序，中间很多公式很烦人的（一定要小心），我写了两天，终于写对了。中间：例如r*(u(i-1,m1,k)+u(i+1,m1,k))形式写成分形式：r*u(i-1,m1,k)+r*u(i+1,m1,k)后面会出错，我也不是很清楚为什么，可能由于舍入误差，或者大数吃掉小数的影响。下面有三个程序具体理论看书，先仔细看书（孙志忠：偏微分方程数值解法）或者网上搜一些理论。Matlab程序：1.function[uu0pexyt]=ADI1(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式，截断)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;fork=1:n+1fori=1:m2+1forj=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfori=1:m2+1forj=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfork=1:n+1fori=1:m2+1u(i,[1m1+1],k)=[exp(0.5*y(i)-t(k))exp(0.5*(1+y(i))-t(k))];u0(i,[1m1+1],k)=[exp(0.5*y(i)-t0(k))exp(0.5*(1+y(i))-t0(k))];endendfork=1:n+1forj=1:m1+1u([1m2+1],j,k)=[exp(0.5*x(j)-t(k))exp(0.5*(1+x(j))-t(k))];u0([1m2+1],j,k)=[exp(0.5*x(j)-t0(k))exp(0.5*(1+x(j))-t0(k))];endendr=h2/(h1*h1);r1=2*(1-r);r2=2*(1+r);fork=1:n%外循环，先固定每一时间层，每一时间层上解一线性方程组%fori=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r2*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+...h2*f(i,2,k);forl=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);forl=1:m1-2%开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1);%回代过程%forl=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendforj=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r2*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);forl=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+...h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);forl=1:m2-2%开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1);%回代过程%forl=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfork=1:n+1fori=1:m2+1forj=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k));%p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k));%e为误差endendend2.function[upexyt]=ADI2(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(D'Yakonov交替方向隐格式)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u'(i,j,k)(引入的过渡层),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;%定义t0是为了f(x,y,t)~=0的情况%fork=1:n+1fori=1:m2+1forj=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));%编程时-t0(k)写成了+t0(k),导致错误;endendend%初始条件fori=1:m2+1forj=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endend%边界条件fork=1:n+1fori=1:m2+1u(i,[1m1+1],k)=[exp(0.5*y(i)-t(k))exp(0.5*(1+y(i))-t(k))];endendr=h2/(h1*h1);r4=1+r;r5=r/2;fork=1:nfori=2:m2u0(i,[1m1+1],k)=r4*u(i,[1m1+1],k+1)-r5*(u(i-1,[1m1+1],...k+1)+u(i+1,[1m1+1],k+1));endendfork=1:n+1forj=1:m1+1u([1m2+1],j,k)=[exp(0.5*x(j)-t(k))exp(0.5*(1+x(j))-t(k))];endendr1=r-r*r;r2=2*(r-1)*(r-1);r3=r*r/2;fork=1:n%外循环，先固定每一时间层，每一时间层上解一线性方程组%fori=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=2*r4*ones(1,m1-1);d(1)=r*u0(i,1,k)+r1*(u(i-1,2,k)+u(i,1,k)+u(i+1,2,k)+...u(i,3,k))+r2*u(i,2,k)+r3*(u(i-1,1,k)+...u(i+1,1,k)+u(i-1,3,k)+u(i+1,3,k))+2*h2*f(i,2,k);forl=2:m1-2d(l)=r1*(u(i-1,l+1,k)+u(i,l,k)+u(i+1,l+1,k)+...u(i,l+2,k))+r2*u(i,l+1,k)+r3*(u(i-1,l,k)+...u(i+1,l,k)+u(i-1,l+2,k)+u(i+1,l+2,k))+2*h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r1*(u(i-1,m1,k)+u(i,m1-1,k)+...u(i+1,m1,k)+u(i,m1+1,k))+r2*u(i,m1,k)+...r3*(u(i-1,m1-1,k)+...u(i+1,m1-1,k)+u(i-1,m1+1,k)+u(i+1,m1+1,k))+2*h2*f(i,m1,k);forl=1:m1-2%开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);end%回代过程%u0(i,m1,k)=d(m1-1)/b(m1-1);forl=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendforj=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=2*r4*ones(1,m2-1);d(1)=r*u(1,j,k+1)+2*u0(2,j,k);forl=2:m2-2d(l)=2*u0(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=2*u0(m2,j,k)+r*u(m2+1,j,k+1);forl=1:m2-2%开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1);%回代过程%forl=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfork=1:n+1fori=1:m2+1forj=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k));%p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k));%e为误差endendend3.function[uu0pexyt]=ADI5(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式，未截断)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;fork=1:n+1fori=1:m2+1forj=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfori=1:m2+1forj=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfork=1:n+1fori=1:m2+1u(i,[1m1+1],k)=[exp(0.5*y(i)-t(k))exp(0.5*(1+y(i))-t(k))];u1(i,[1m1+1],k)=[exp(0.5*y(i)-t0(k))exp(0.5*(1+y(i))-t0(k))];endendr=h2/(h1*h1);r1=2*(1-r);r2=r/4;r3=2*(1+r);fork=1:nfori=2:m2u0(i,[1m1+1],k)=u1(i,[1m1+1],k)-r2*(u(i-1,[1m1+1],k+1)-...2*u(i,[1m1+1],k+1)+u(i+1,[1m1+1],k+1)-u(i-1,[1m1+1],k)+...2*u(i,[1m1+1],k)-u(i+1,[1m1+1],k));endendfork=1:n+1forj=1:m1+1u([1m2+1],j,k)=[exp(0.5*x(j)-t(k))exp(0.5*(1+x(j))-t(k))];endendfork=1:n%外循环，先固定每一时间层，每一时间层上解一线性方程组%fori=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r3*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+...h2*f(i,2,k);forl=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);forl=1:m1-2%开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1);%回代过程%forl=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendforj=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r3*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);forl=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+...h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);forl=1:m2-2%开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1);%回代过程%forl=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfork=1:n+1fori=1:m2+1forj=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k));%p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k));%e为误差endendend[upexyt]=ADI2(0.01,0.001,100,100,1000);surf(x,y,e(:,:,1001))t=1的误差曲面下面是三种 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 的误差比较：1.[uu0pexyt]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2）代替）（t=1时的误差）2.[uu0pexyt]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)3.[upexyt]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式)surf(x,y,e(:,:,11))( 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示t=1时的误差）下面是相关数据：1:[uu0pexyt]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2)代替）e(:,:,11)=Columns1through600000000.000409470.000251820.000190770.000171120.0001760400.000573590.000429710.000354020.000325650.0003362800.000662360.000546890.000474080.000445960.0004626700.000721520.000620010.000550810.000524420.0005455300.000761640.00065760.000585220.000557320.0005798400.000783360.000659930.000575570.000541610.0005620900.000781610.000618720.000516460.000474290.0004896400.000736210.00051480.000399790.000354390.0003631300.000569640.000316880.000220510.00018840.000191920000002.[uu0pexyt]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)e(:,:,11)=Columns1through600000000.000270060.000163050.000121040.00010710.0001099500.000377540.000278170.00022530.000204830.0002111600.000435390.000353860.000302070.000281240.000291400.000473980.000401040.000351130.000331110.0003440500.00050030.000425350.000373090.00035190.0003657100.000514790.000426990.000366810.000341640.000354100.000514150.000400560.000328870.00029850.0003076400.000485040.00033350.000254110.00022210.0002270600.000376090.000205320.000139560.000117180.000119020000003.[upexyt]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式)e(:,:,11)=Columns1through600000008.6469e-0061.4412e-0051.8364e-0052.091e-0052.2174e-00501.4412e-0052.4777e-0053.2047e-0053.6716e-0053.8961e-00501.8364e-0053.2047e-0054.1789e-0054.8054e-0055.1008e-00502.091e-0053.6716e-0054.8054e-0055.5353e-0055.8764e-00502.2174e-0053.8961e-0055.1008e-0055.8764e-0056.2389e-00502.2118e-0053.8698e-0055.0523e-0055.8126e-0056.171e-00502.055e-0053.5581e-0054.6157e-0055.2942e-0055.6197e-00501.707e-0052.8951e-0053.7128e-0054.2365e-0054.4952e-00501.0851e-0051.7698e-0052.2265e-0052.5203e-0052.672e-0050000001.[uu0pexyt]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2)代替）Columns7through11000000.000203480.000262280.000383380.0006600800.000386070.000483210.000647170.0009166800.000526350.000642030.000816370.001051700.00061740.000742720.000921110.001141700.000656510.000789640.000977240.001205100.000640510.000781160.000985940.001243300.000564740.000708220.000933320.001247800.000425470.000555260.000786160.001184400.000227350.000309460.000490040.000924020000002.[uu0pexyt]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)Columns7through11000000.000128260.000167980.000249860.0004363700.000244440.000310230.000421730.0006051300.000334010.000412570.000531790.0006935800.000392160.000477420.000599860.0007526300.000417040.000507610.000636420.0007943900.000406570.00050210.000642260.0008198400.000357840.0004550.000608280.0008233400.000268660.000356280.000512630.000782400.000142620.000197890.000319560.00061130000003.[upexyt]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式)Columns7through11000002.2118e-0052.055e-0051.707e-0051.0851e-00503.8698e-0053.5581e-0052.8951e-0051.7698e-00505.0523e-0054.6157e-0053.7128e-0052.2265e-00505.8126e-0055.2942e-0054.2365e-0052.5203e-00506.171e-0055.6197e-0054.4952e-0052.672e-00506.1116e-0055.5803e-0054.4817e-0052.6785e-00505.5803e-0055.1239e-0054.1529e-0052.5153e-00504.4817e-0054.1529e-0053.42e-0052.126e-00502.6785e-0052.5153e-0052.126e-0051.3869e-005000000友情提示： 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 范本是经验性极强的领域，本范文无法思考和涵盖全面，供参考!最好找专业人士起草或审核后使用。