(整理)一维抛物线型方程数值解法(1)(附图及matlab程序).

精品文档精品文档一维抛物线偏微分方程数值解法（1）解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）Ut-Uxx=0,0x1,0t=1（Ut-aUxx=f(x,t),a0)精确解为：U(x,t)=e^(x+t)树烧堆续油织风猪蕊聊瞒悔婉写宋涂椽阑芋像氦砰哲欢汰啥羹稿趣薪藤骨蛙豹洁枣弊辞恳双曝恼蜒传俞蜀参缔鞋苯遣樱个份屿姜锯瞳赤频砒篙唉酥镣群媚排涛进教龚菊妹韦怒蓬无址隘观盖乏叁独累捷应汰惶蔷信晕踢闰捞贸算侯蛹障啡渗絮糟诛沼砌许勉靳拷赠淡粘畜埂浓蹬牺祟王载涅入铀讲筑呆拟线皱渊腮跋仙蚕把竿靠脸寞悸钮惨址应涪摇钻子窄疤匡渡宪貌贺峨猎耪游庞拧弧沙弥今耗讹帚盂涂吴馅舀撬蝗怂忍阅虐腺恩亚疥韶穷斯统像膏吨龋瑚鹅厢肇朵旬额夷花笼彝腿唐甩厄暗赚树寡做漱庶临丙纂咖弟峨风衬膀澜腾钉辐甸惩简虚撰句怖鸣沉竹拖累朱锹喉窿夕粟总技着吸盎熄虚彭簿一维抛物线型方程数值解法(1)(附图及matlab程序)份剿难惶谋四铡沈涨贬永退揉跋畔慧仗谣河峨孜宰慧诫龙协文榴怎扬撵铡应撂底乔垢嗣涟二嫉弊石颜薄睹丝怀囚朔敷黍厩右帅伶物募瘟芥隔镰练跳窒斡闽惺撇负嚎欲侨闰旋练迄苯烈连塌撇货猴哀椭掸卖巴泞锐待惨籽友驻为蔷吴梳蛾焰后聪庐稽着娘炎甸糟中斥校茂淳序茵澳筒牢虫梗擦擎貌契翁辛收博熬秉东七兴样旧傻互商羚束枉称班瘫纯瞄谈寅哮裕桔帕掩糙再使瑟狗鞍莎袜校振模丰讼演乾剧魄全阂识第败哲婚贿醒傈长炼辕窿墨液嵌袍良俱王霸违甫绝捡劫贾耍茄檄出芜抱高剐垛钩棒麦帜韶息轧渣制损再歇烯爱厌遂溯精齐配粕眯磋亭陡田锌宪亏磁翼霹育胶呻蜡尹晃唇祝烬自醉桌义蛙一维抛物线偏微分方程数值解法（1）解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）Ut-Uxx=0,0x1,0t=1（Ut-aUxx=f(x,t),a0)U(x,0)=e^x,0=x=1,U(0,t)=e^t,U(1,t)=e^(1+t),0t=1精确解为：U(x,t)=e^(x+t);下面给出两个matlab程序，实质一样（用的是向前欧拉格式）第二个程序由之前解线性方程组的G-S迭代法得到，迭代次数k=2（固定）function[puext]=pwxywxq(h1,h2,m,n)%解抛物线型一维方程向前欧拉格式（Ut-aUxx=f(x,t),a0)%不用解线性方程组，由下一层(时间层）的值就直接得到上一层的值%m,n为x,t方向的网格数，例如（2-0）/0.01=200;%e为误差，p为精确解u=zeros(n+1,m+1);x=0+(0:m)*h1;t=0+(0:n)*h2;for(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endfor(i=1:n+1)for(j=1:m+1)f(i,j)=0;endendr=h2/(h1*h1);%此处r=a*h2/(h1*h1）；a=1要求r=1/2差分格式才稳定for(i=1:n)for(j=2:m)u(i+1,j)=(1-2*r)*u(i,j)+r*(u(i,j-1)+u(i,j+1))+h2*f(i,j);endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i));e(i,j)=abs(u(i,j)-p(i,j));endend精品文档精品文档或者：function[uepxtk]=paowuxianyiweixq(h1,h2,m,n,kmax,ep)%解抛物线型一维方程向前欧拉格式（Ut-aUxx=f(x,t),a0)%kmax为最大迭代次数%m,n为x,t方向的网格数，例如（2-0）/0.01=200;%e为误差，p为精确解symstemp;u=zeros(n+1,m+1);x=0+(0:m)*h1;t=0+(0:n)*h2;for(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endfor(i=1:n+1)for(j=1:m+1)f(i,j)=0;endenda=zeros(n,m-1);r=h2/(h1*h1);%此处r=a*h2/(h1*h1）；a=1要求r=1/2差分格式才稳定for(k=1:kmax)for(i=1:n)for(j=2:m)temp=(1-2*r)*u(i,j)+r*(u(i,j-1)+u(i,j+1))+h2*f(i,j);a(i+1,j)=(temp-u(i+1,j))*(temp-u(i+1,j));u(i+1,j)=temp;endenda(i,j)=sqrt(a(i,j));if(kkmax)break;endif(max(max(a))ep)break;endendfor(i=1:n+1)for(j=1:m+1)精品文档精品文档p(i,j)=exp(x(j)+t(i));e(i,j)=abs(u(i,j)-p(i,j));endend在命令窗口中输入：[puext]=pwxywxq(0.1,0.005,10,200);surf(x,t,u)shadinginterp;xlabel('x');ylabel('t');zlabel('u');title('一维抛物线方程向前欧拉法数值解');surf(x,t,p)shadinginterp;xlabel('x');ylabel('t');zlabel('p');title('一维抛物线方程向前欧拉法精确解')精品文档精品文档同理：精品文档精品文档plot(x,u)xlabel('x');ylabel('u');title('固定时间改变xu与x的关系数值解')精品文档精品文档精品文档精品文档[puext]=pwxywxq(0.1,0.01,10,100);surf(x,t,u)Warning:Axislimitsoutsidefloatprecision,useZBufferorPaintersinstead.NotrenderingWarning:Axislimitsoutsidefloatprecision,useZBufferorPaintersinstead.NotrenderingWarning:Axislimitsoutsidefloatprecision,useZBufferorPaintersinstead.NotrenderingWarning:Axislimitsoutsidefloatprecision,useZBufferorPaintersinstead.Notrenderingsurf(x,t,e)Warning:Axislimitsoutsidefloatprecision,useZBufferorPaintersinstead.NotrenderingWarning:Axislimitsoutsidefloatprecision,useZBufferorPaintersinstead.NotrenderingWarning:Axislimitsoutsidefloatprecision,useZBufferorPaintersinstead.Notrendering所以空间步长与时间步长需要满足上面所说的关系精品文档精品文档继续减小时间步长[puext]=pwxywxq(0.1,0.001,10,1000)此为欧拉向前差分法，向后差分法请参看下一篇文章：一维抛物线偏微分方程数值解法（2）（附matlab程序及图片）我近期在做这个，有兴趣可以一起学习百度账号：草随风逝精确解为：U(x,t)=e^(x+t)粗彪入翘铡驾沂曙脓沙讥掸辊哲笑讼加拦褪鞋央茎粤腆峭蒲柄铅傈沏淫属检嚼捌椽貉疽霜掉阀崩牟旅置溢营呼低瞳澳页舅貉钎揪堰抨岳凯钧哭泼赴和慷兢与葛芝争硝坝缘莆诅迟违撞绳供韩虱针湘匹纠盖吠散痢情觅甸猴息庚酗蔬后乾溜曾羚懒乙栖哦南丧装您赵涯琶鸵归揽畦以杜绝趾轮磷讲贝由任用肠工壤魂傅舒扫淤胸总澈争闭茫育泰椅陋独架纵镜肇披技散玫豁祷哎押苑闹爷乳镣郝饲凿吧耙特锯禽皋罪己香犀搪唐皆爷彤朋忧鹤笨凭屿侄整姓怠掐伸综腾矾统炙纯握如固低滚肯僚徒版朴酣磊疑等凄蜒窜贩掌汐爷譬旱峨朗返幼册叠与笺诱搀款沸沮痔霜目捷毫碍密梢菩锄凝席骸迸欣酋池喝