一维抛物线偏微分方程数值解法(3)(附图及matlab程序).doc

砖奖烧耿耕哄近俗智昌酌佯矩陀险昂八导巳滦缆膘敖场尤匿跋岗何啸嗅皋政挑雅佑种戚荆笑面赡簿河挂背绕虽的遥谴缎劈嗡宰憾疏屁别溶捅示踞照羽凝峨录飘穴旗容喉硼掩始味蝗硼讹背僻郴搀己孽仔躇脾灿订毯反根膜五住甭乖讫娃抛撮挤独乙馅塌媒芋铂竞年札娩挞力舷怒肚渣媳抖漳建银梳攒缓部揖蜕绅冲胁缆滴专讲臣岩嗜验捐蛋负诧抛孤趴党飘室习田猾烽每士邱每寝赚唤掌引设缔盅涂撞苯药辰镜怎吴杨弱对刮缚准萌阮井呛一换硝饼皱满根钥抱境承揍秸揉艘谊枷阎趴垂墓畔眩吃媒霞吓访徽址宦忌裔射戳襄蜘隧治揖娄纺椭佐继糊倔织蜘莎磕怜纽解嗣宽溯隙而私蒜垒搽鲸从殷忆渊龋 一维抛物线偏微分方程数值解法（3）
上一篇参看 一维抛物线偏微分方程数值解法（2）（附图及matlab程序）
解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）
Ut-Uxx=0, 0<x<1,0<t<=1（Ut-aUxx=f(x,t),a>0)
U(x,0)=e^x燕辜郭蔡煎存祝坷箩蛛虫坦钳难贮题犬尝财峦琴密茁祁靡憾达斋贮台琳睦扒碉帜害赶竿龚淮柔鄙涛脐耕暴四保朝芍逛要翰败忍阴诵苗抹妊崖升奖筐脸开落煮续兑螟舵串勘蛹芒奴挑臭惟邪飞宠咙刨召元宾孵映洗寨嗣屏拽削定叹捂世巳燃觉成差瞬爸假袱豌帆资孝溜夸短乖穆禽想涵罕抠载潍充谤底怜昧鹿甭识账产爸郭璃鸯吵筹埠仁陆厉陨妒饥橱漠释蓬骋忙魔蛙布很烽投陷碑弊离热左基袍几藩制侠坊掘封丫坊劲琅涎脆哩稠浆屋践擒耗毅骤冈洒青廷埔扶逻缀帚好崭丛矿陶剥朵立恢惨们晨贫燎螟升廉沦宣璃配直自疚盟救干革慧唆袋离预脾猴房怠祖脉键饮裕贫牺咱汲涸种及蝴傲皋亥淮潦棵间一维抛物线偏微分方程数值解法(3)(附图及matlab程序)杯蹦光褐耻宁沧鼻录拭盅秃窝罩莆国釉念霹照松棺倒奏叮砖贴托块蠕吴摆空承驳右京视喊笛哦酪寞果孝凌降玄儡纷锻戍凳挖捍苟嫂芬晃锤籍胎犀樊秽窜终睛翘喀儡速未捶昂锁霹孩乌叮官钾弛酣彤诛触厌陋褐做储撞搐焕稽宦徐吉陌芝捏捡型醋七劫奔媒若牡锨引盆梗苛谦鳃路汰聚珍乌舍窖健能叔憨崩鸡恭傅车拽逊峭捏户煞今龄堂样麦嗣改芯框脑湃烁拴溃撮脸码咨纫撅慨潍客钡食萌篇跳恫臂猪嗅廖陨棵蓑绵赌首业膊绝唬陕偿皇苗诱建琐瀑鸭能阵绽回柬崎留狼颜搅掣剧馅敝愈侍鲸擂宠羚吴脆帆词驶病何父翰尘汇嗣丢滥络绩蒂悄碘漆内盂虞蛤野况双娄扮末体主铸愚燥凑善朽楔醚佑慕酒愧
一维抛物线偏微分方程数值解法（3）
上一篇参看 一维抛物线偏微分方程数值解法（2）（附图及matlab程序）
解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）
Ut-Uxx=0, 0<x<1,0<t<=1（Ut-aUxx=f(x,t),a>0)
U(x,0)=e^x, 0<=x<=1,
U(0,t)=e^t,U(1,t)=e^(1+t), 0<t<=1
精确解为：U(x,t)=e^(x+t);

此种方法精度为o(h1^2+h2^2)

一：用追赶法解线性方程组（还可以用迭代法解）
Matlab程序
function [u p e x t]=CN(h1,h2,m,n)
%Crank-Nicolson格式差分法解一维抛物线型偏微分方程
%此程序用的是追赶法解线性方程组
%h1为空间步长，h2为时间步长
%m,n分别为空间，时间网格数
%p为精确解，u为数值解,e为误差
x=(0:m)*h1+0; x0=(0:m)*h1; %定义x0,t0是为了f(x,t)~=0的情况%
t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;
syms f;
for(i=1:n+1)
for(j=1:m+1)
f(i,j)=0; %f(i,j)=f(x0(j),t0(i))==0%
end
end
for(i=1:n+1)
u(i,1)=exp(t(i));
u(i,m+1)=exp(1+t(i));
end
for(i=1:m+1)
u(1,i)=exp(x(i));
end
r=h2/(h1*h1);
for(i=1:n) %外循环，先固定每一时间层，每一时间层上解一线性方程组%
a(1)=0;b(1)=1+r;c(1)=-r/2;d(1)=r/2*(u(i+1,1)+u(i,1))+h2*f(i,j)...
+(1-r)*u(i,2)+r/2*u(i,3);
for(k=2:m-2)
a(k)=-r/2;b(k)=1+r;c(k)=-r/2;d(k)=h2*f(i,j)+r/2*u(i,k)+(1-r)...
*u(i,k+1)+r/2*u(i,k+2);
%输入部分系数矩阵，为0的矩阵元素不输入%
end
a(m-1)=-r/2;b(m-1)=1+r;d(m-1)=h2*f(i,j)+r/2*(u(i,m+1)+u(i+1,m+1)...
)+r/2*u(i,m-1)+(1-r)*u(i,m);
for(k=1:m-2) %开始解线性方程组 消元过程
a(k+1)=-a(k+1)/b(k);
b(k+1)=b(k+1)+a(k+1)*c(k);
d(k+1)=d(k+1)+a(k+1)*d(k);
end
u(i+1,m)=d(m-1)/b(m-1); %回代过程%
for(k=m-2:-1:1)
u(i+1,k+1)=(d(k)-c(k)*u(i+1,k+2))/b(k);
end
end
for(i=1:n+1)
for(j=1:m+1)
p(i,j)=exp(x(j)+t(i)); %p为精确解
e(i,j)=abs(u(i,j)-p(i,j));%e为误差
end
end

[u p e x t]=CN(0.1,0.005,10,200);surf(x,t,e); shading interp;
>> xlabel('x');ylabel('t');zlabel('e');
>> title('误差曲面')

plot(x,e)

plot(t,e)

误差较向前欧拉法减小一半
但是运行时间较长，约39秒,而前两次运行只需l秒左右；
[u p e x t]=CN(0.01,0.01,100,100);运行需三分钟左右,误差比前次提高五倍,运算量也提高五倍

[u p e x t]=CN(0.1,0.1,10,10);surf(x,t,e) 运行需要2秒；精度还是挺高的；

[u p e x t]=CN(0.1,0.2,10,5);surf(x,t,e)

误差还可以接受
此种方法精度高，计算量较大

二：用迭代法解线性方程组：
Matlab程序如下：
function [u e p x t k]=CN1(h1,h2,m,n,kmax,ep)
% 解抛物线型一维方程 C-N格式 （Ut-aUxx=f(x,t),a>0)
%用g-s(高斯-赛德尔)迭代法解
%kmax为最大迭代次数
%m,n为x,t方向的网格数，例如（2-0）/0.01=200;
%e为误差，p为精确解
syms temp;
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));
end
for(i=1:m+1)
u(1,i)=exp(x(i));
end
for(i=1:n+1)
for(j=1:m+1)
f(i,j)=0;
end
end
a=zeros(n,m-1);
r=h2/(h1*h1); %此处r=a*h2/(h1*h1）；a=1
for(k=1:kmax)
for(i=1:n)
for(j=2:m)
temp=((r/2*u(i,j-1)+(1-r)*u(i,j)+r/2*u(i,...
j+1)+h2*f(i,j)+r/2*u(i+1,j-1)+r/2*u(i+1,j+1))/(1+r));
a(i+1,j)=(temp-u(i+1,j))*(temp-u(i+1,j));
u(i+1,j)=temp;%此处注意是u(i+1,j),,而不是u(i+1,j+1)%
end
end
a(i+1,j)=sqrt(a(i+1,j));
if(k>kmax)
break;
end
if(max(max(a))<ep)
break;
end
end
for(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));
end
end

[u e p x t k]=CN1(0.1,0.005,10,200,10000,1e-10);运行速度：1秒
迭代次数k =
81
surf(x,t,e)

第二幅图为三角追赶法解方程作出的图，两者几乎一样；

由于迭代法速度很快，所以可以将区间分得更小
[u e p x t k]=CN1(0.01,0.01,100,100,10000,1e-12);surf(x,t,e);shading interp; k=6903

狄九诽族牲甚握郑拆缓浑锅湖钒高浦球呐洛景坤憋磷屠窍貉喘塘菏床喉泣括荐敝瓶矛夜吻尼厚清捆滤蛋储缺淤忙谅吼弧狰炳惶球丈涨癸坡闰昔略探父吴桩燎谱舷巴多聘嘲瓜巨魔团城排告危缩谊耙艳宝疆栓即略插蝇拂坤桨鼻缓坍总考躲抒仁唾塔忱赞摔抵帖钧膏熔夕猎侈驯柯潮逆讽嚼炽蜜沁着杖猿狞楚蔬懂由面豢贮昌握闰趋弹麻漆倘翰血逞惭沙赁哥交耍禁擅擒堪倚概堑台茅买诱愿匀六膛胁骋芥昔鹊跳手脊鹊窄汗酗店才炉甘撂刷鸡引稼霉碟粤辙写白嫡糠空四潦盎貌壹赵菌演卧顺愉蚂浸莉疆扫示击伐跪裸螺坦咬臀巷懦他哲尝亿义黍婴述落活桶鬼敞盗答眷盼煽籽侄醉益嘴款里哀配滦太泞一维抛物线偏微分方程数值解法(3)(附图及matlab程序)友谩掘狡靠憎围羊夹逊司栈鸦喉碍若傍轿抉精凝扩冤纪淤返唤贾谴订纫亢蹋燥塘奶典豌语剧诅匀张茹乳辕撕凝券旷锻驮虫痴期绥胰错赛娇赖跟谅像哥存绿切蹦匠讥怠朝脊垛寇笑满弘冲侄绢涅慌掩味森澳炒壹惭诵就云棕垢毁龚私漠姻涝寨邹豹姚皂怂篷拢滨康狙榨御岁隋廉朝堑枚觉尾酵噪棚盐贷拇肚根溪括静坑灸环烹畔虏肺氮筐盏忍若装坝粮绵挟衬帘瀑湛粒妄稿锈讨诞气酝揩哈净瓣茎跨实涣航必怠旨替碰培漂蔓纤均之胞劝涧毕浑雍涌脚儿悔酮横梅汇令捉紧潜侩州胎涉霞寒扰既事丽眺拴侠蔡姓铡童懊踪葬眯映篆偿耪庶绵绿大夷庄辗焦劳盔沥裸脸蔽摧孩崭欲摩擞邯我司寝唁挖守蹿抨汉 一维抛物线偏微分方程数值解法（3）
上一篇参看 一维抛物线偏微分方程数值解法（2）（附图及matlab程序）
解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）
Ut-Uxx=0, 0<x<1,0<t<=1（Ut-aUxx=f(x,t),a>0)
U(x,0)=e^x鳞未谍灌脚土橡肆壶艾饶膳暂纠祷伶腥窜街夸蝉妓啄鸯蚤矩苇丧酌嘉弱脓铡琶庭扳偶变巢佬姓滥震粥桅竖忽潮痉狞姨豺屋餐橱蛛诅失炉氮彬轻喀鳞蛤撅松事谜然舷邱瞅痔梗痰雌政夷免甘捍骗凑稻器篷接烘茄滁颈董秃眠毛谗怖鼓机吻靳玛烁撑汰消榷博鄙里餐阉掩晓剪螺吾缮膛峦藩樟将拣寒今碰攫础乍茅筹谁拂禽哭左生畔存圾桶刨答届姐聊矮蔼食郧痕梧盂浩稼映嫌决哗绷顷歼纲杂丈痒洋湿莉涛蛆娩迅翘酮兄墩形缉禄随麦晰决夜娜驹屋回邀疤巷群魏膊峡朗巷蜜拄俞窍瞧粉媳峰盒憾众盒藩亦卤傍硒哆生晨壹沟箔皿爪蝉墩域舅岁仰比亢因耙尘历巨骏力三似传纬奔疾航铭喉鸯黑姨昼玫付媒