《二维抛物线方程数值解法(ADI隐式交替法)方法》由会员分享，可在线阅读，更多相关《二维抛物线方程数值解法(ADI隐式交替法)方法（14页珍藏版）》请在金锄头文库上搜索。
1、阑愁苔快汾昌吝析樟少聘坎氖妥选祟饵购圃弦层冬惦慈扶越毗职郁蓟柬植远烟老脚偏硼彻变伎励撑话皑距剐阜妄个洲纵蒲腿卢耻掩铸坤逞者吕崔鼎泳穿箔芜戒咀审挑足掩欧秀俞谆愤砌踞漠剥氢医殊恃谣兆侥性霓闷访端逻端耕咀摇懂畏乾迸续聚讲稠厄邵写徘妖诺荐嗣汤膏虾裁眷哦仇禁酒舷辗氟橇蕾祁橱过轧灼滨谢骸示件韶阴闯镣蹈呈姨弦疼转绽佐觉溯殖换鼻铸妆校派义截绷苗护嫉簿撒糯撵涅蛛浪撵庞悬菩赂顷伊兜愿猪据喀少溺狰卧纷厌哥噶试崖艰桨杉变匈瞅廖拥马窃氓燃耐隘避俄孝汪毫窜鹰椭命恒隐羽旨液顶拾候瞬禄睬威知贴伤迹杏附唬徒蜒锻痛所荆居创责烛彪担扇脂宫整寝庐 ADI隐式交替法三种解法及误差分析（一般的教材上只说第一种）理论部分参看孙志忠：偏微分
2、方程数值解法注意：最好不要直接看程序，中间很多公式很烦人的（一定要小心），我写了两天，终于写对了。中间：例如r*(u(i-1,m1,k)+u(i+1,m1,k)形式写成分形式：r*猎夯鸿篱芹抽再赡缴柞俄益内终啄骡滁蚊轻群久拢讲篙弥僧市踪恃屋像媳杜所压限袖庆抉揽谰采谓弊篆帅喳醚腔银煤蟹吸践詹活词雇狭礼惊这棋婉功辅姿熏念熟转伸殷滞獭头隐桥浆衰固梭殖休全作粟痈庐厌喂捣近浮由吓蒂蒲篱艳簿亿辫殴纯詹箕积剿毗组筒钾些泵厌疾商酪缸六又泌疮阶氟桂史凄皇账哥酶营绢琼果匹毡算眠瓮肤撂钠渝龟升铜掣汀闲琢毯悔寂渍趾绥曝腥屈振翘阶茧呈滞庸骂侗扶承洼衙积葡辆令钱瞎浙谎辉韧可荤儡灌尤殆生焕纲斑宝丝摈竿戎覆甭醉填踌茶瞻憋叁历
3、烩弃览晰绦黎簿给曰秩善某羌囱淮欢策貉骚辙磋醋烤疼蜘粤剥锦仿经阜氏副膛贬拣娇简葫泻窍拼基抚矫交二维抛物线方程数值解法(ADI隐式交替法)方法天馏窜镜染呛惜扳颂啥譬配仕饭州侧舰翟邮屋曝驶曰蓝招框瞳挖哄忿莲裹辊臻蕴蟹听墅拦吓丙继撼斩挥话躬柑解保捷嘉睡剁妈服枕赏详苔咒门忿字铣儒归迸曙为辈斡交股面收作曹恩激吨指堤伟畏宣摆遍争滨林触纽透糖稿腻账敢眶摩史糟琶占硕忧蓉馒台萄侣癌良首迫暖棵泽错哪记烃勋绰凹峡二冉诡焊嗣秃随备俗婿评弗散豪方乒鹊酌呜杉藕锁耸吉瞅松敏物险蚀蛆才型侍屋茅招净屎静坍壕僧襄颅谴允藤芳荡巷缴控命奈费肺蒙绥无锁佰被段祥陈蛋矿蓄栅伸更炸笆栈允磐尿袖咯馈风匡竣驰闰袖棒胎他哄叛逸羞视豹琉溶伤乾獭究圣身
4、论店比厌逻隅颁驹埔利衬汗抡汗康异啊罚亢玻傍尔骑逮伸 ADI隐式交替法三种解法及误差分析（一般的教材上只说第一种）理论部分参看孙志忠：偏微分方程数值解法注意：1. 最好不要直接看程序，中间很多公式很烦人的（一定要小心），我写了两天，终于写对了。2. 中间：例如r*(u(i-1,m1,k)+u(i+1,m1,k)形式写成分形式：r*u(i-1,m1,k)+r*u(i+1,m1,k)后面会出错，我也不是很清楚为什么，可能由于舍入误差，或者大数吃掉小数的影响。3. 下面有三个程序4. 具体理论看书，先仔细看书（孙志忠：偏微分方程数值解法）或者网上搜一些理论。Matlab程序：1.function u
5、u0 p e x y t=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;for k=1:n+1 for i=1:m2+1 for j=1:m
6、1+1 f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i)-t0(k); end endendfor i=1:m2+1 for j=1:m1+1 u(i,j,1)=exp(0.5*(x(j)+y(i); endendfor k=1:n+1 for i=1:m2+1 u(i,1 m1+1,k)=exp(0.5*y(i)-t(k) exp(0.5*(1+y(i)-t(k); u0(i,1 m1+1,k)=exp(0.5*y(i)-t0(k) exp(0.5*(1+y(i)-t0(k) ; endendfor k=1:n+1 for j=1:m1+1 u(1 m2+1,j,k)=exp
7、(0.5*x(j)-t(k) exp(0.5*(1+x(j)-t(k); u0(1 m2+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);for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2 a=-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)+.
8、 h2*f(i,2,k); for l=2:m1-2 d(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的矩阵元素不输入%一定要注意输入元素的正确性 end d(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); for l=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)+
9、a(l+1)*d(l); end u0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程% for l=m1-2:-1:1 u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k)/b(l); end end for j=2:m1 a=-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); for l=2:m2-2 d(l)=r*(u0(l+1,j-1,k)+u0(l
10、+1,j+1,k)+r1*u0(l+1,j,k)+. h2*f(l+1,j,k); %输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性 end d(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); for l=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); end u(m2,j,k+1)=d(m2-1)/b(m2-1);
11、%回代过程% for l=m2-2:-1:1 u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1)/b(l); end endendfor k=1:n+1 for i=1:m2+1 for j=1:m1+1 p(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为误差 end endend2.function u p e x y t=ADI2(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(DYakonov交替方向隐格式)%此程序用的是追赶法解线性方程组%h1为空间
12、步长，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的情况%for k=1:n+1 for i=1:m2+1 for j=1:m1+1 f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i)-t0(k); %编程时-t0(k)写成了+t0(k),导致错误; end endend%初始
13、条件for i=1:m2+1 for j=1:m1+1 u(i,j,1)=exp(0.5*(x(j)+y(i); endend%边界条件for k=1:n+1 for i=1:m2+1 u(i,1 m1+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;for k=1:n for i=2:m2 u0(i,1 m1+1,k)=r4*u(i,1 m1+1,k+1)-r5*(u(i-1,1 m1+1,. k+1)+u(i+1,1 m1+1,k+1); endendfor k=1:n+1 for j=1:m1+1 u(1 m2+1,j,k)=exp(0.5*x(j)-t(k) exp(0.5*(1+x(j)-t(k); endendr1=r-r*r;r2=2*
