拋物型微分方程

數(shù)學(xué)與計(jì)算科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)工程名稱拋物型方程數(shù)值解所屬課程名稱微分方程數(shù)值解法實(shí)驗(yàn)類型驗(yàn)證實(shí)驗(yàn)日期班級(jí)信計(jì)0902學(xué)號(hào)姓名成績一、實(shí)驗(yàn)概述：【實(shí)驗(yàn)?zāi)康摹?熟練掌握向前差分格式和Crank-Nicolson格式2求解系數(shù)矩陣為三對(duì)角陣的方程組的追趕法【實(shí)驗(yàn)原理】【實(shí)驗(yàn)環(huán)境】MicrosoftVisualC++6.0二、實(shí)驗(yàn)內(nèi)容：【實(shí)驗(yàn)方案】用向前差分格式和Crank-Nicolson格式求解一維熱傳導(dǎo)方程初邊值問題0<x<1,0<t<Tu(x,0)=sin(pi*x)0<x<1u(0,t)=u(1,t)=00<t<T【實(shí)驗(yàn)過程】〔實(shí)驗(yàn)步驟、記錄、數(shù)據(jù)、分析〕1，局部算法代碼：(1)三角分解voidfenjie(doublea,doubleb,doublec,doublee[m-1],doubled[m-1])//三角分解{ inti;d[0]=a; for(i=2;i<=m-1;i++) { e[i-1]=b/d[i-2]; d[i-1]=a-c*e[i-1]; }}voidhuidai(doubleg[m-1],doubley[m-1],doublee[m-1],doubled[m-1],doublec,doublex[m-1])//回代求解{ inti; for(i=2,y[0]=g[0];i<=m-1;i++) y[i-1]=g[i-1]-e[i-1]*y[i-2];//追的過程 for(i=m-2,x[m-2]=y[m-2]/d[m-2];i>=1;i--) { x[i-1]=(y[i-1]-x[i]*c)/d[i-1];//趕的過程 }}(2)Crank_Nicolson法voidCrank_Nicolson(doubler,doubleU[m+1][n+1]){doublea,b,c,g[m-1],e[m-1],d[m-1],y[m-1];doublexx[m-1];inti,j;a=1+r;b=-1.0/2*r;c=b; for(j=0;j<=n;j++) {for(i=0;i<m-1;i++) { if(i==0) g[0]=(1-r)*U[1][j]-b*U[2][j]-b*(U[0][j]+U[0][j+1]); if(i==m-2) g[m-2]=(1-r)*U[m-1][j]-b*U[m-2][j]-b*(U[m][j]+U[m][j+1]); else g[i]=(1-r)*U[i+1][j]-b*(U[i+2][j]+U[i][j]); } fenjie(a,b,c,e,d);huidai(g,y,e,d,c,xx); for(i=0;i<m-1;i++) U[i+1][j+1]=xx[i]; } }(3)目標(biāo)函數(shù)及解析解doubleu(doublet,doublex)//微分方程的精確解{ doubleuu; uu=exp(pow(pi,2)*(-t))*sin(pi*x); returnuu;}2，結(jié)果分析：根據(jù)上面的實(shí)驗(yàn)原理，可以得出：在每一時(shí)間層，需要求解的隱式差分方程形成了一個(gè)線性代數(shù)方程組，它的系數(shù)矩陣是三對(duì)角形矩陣，即僅在主對(duì)角線及其相鄰二條對(duì)角線上有非零元素，故可以使用三角分解法。使用Crank_Nicolson隱式法，每時(shí)間層包含較多的計(jì)算工作量，但是其優(yōu)點(diǎn)為其穩(wěn)定性要求對(duì)步長比的限制大為放寬，而這正是我們所期望的。【實(shí)驗(yàn)結(jié)論】〔結(jié)果〕下面僅顯示X軸上取x=0.5處時(shí)間軸T的一列間斷點(diǎn)出的函數(shù)值：【實(shí)驗(yàn)小結(jié)】〔收獲體會(huì)〕通過本次實(shí)驗(yàn)，我重新復(fù)習(xí)了三角分解法，也對(duì)課本上的Crank_Nicolson方法有了更深的了解，同時(shí)也提高了自己的上機(jī)編程能力，感受到了理論與實(shí)踐相結(jié)合的樂趣。三、指導(dǎo)教師評(píng)語及成績：評(píng)語評(píng)語等級(jí)優(yōu)良中及格不及格1.實(shí)驗(yàn)報(bào)告按時(shí)完成,字跡清楚,文字表達(dá)流暢,邏輯性強(qiáng)2.實(shí)驗(yàn)方案設(shè)計(jì)合理3.實(shí)驗(yàn)過程〔實(shí)驗(yàn)步驟詳細(xì),記錄完整,數(shù)據(jù)合理,分析透徹〕4實(shí)驗(yàn)結(jié)論正確.成績：指導(dǎo)教師簽名：批閱日期：附錄1：源程序#include<stdio.h>#include<stdlib.h>#include<math.h>#definen800#definem20voidfenjie(doublea,doubleb,doublec,doublee[m-1],doubled[m-1]);//三角分解voidhuidai(doubleg[m-1],doubley[m-1],doublee[m-1],doubled[m-1],doublec,doublex[m-1]);//回代求解voidCrank_Nicolson(doubler,doubleU[m+1][n+1]);//追趕法求解方程組doubleu(doublet,doublex);//精確解voidmain(){doublex[m+1],t[n+1]; doubleU[m+1][n+1]; doubler,h,k;//步長 inti,j;//賦初值 h=1.0/m; r=1/sqrt(20); k=pow(h,2)*r; for(i=0;i<=m;i++) x[i]=i*h; for(j=0;j<=n;j++)t[j]=j*k; for(i=0;i<=m;i++)//初值 { U[i][0]=sin(pi*x[i]); } for(j=0;j<=n;j++)//邊值 { U[0][j]=0; U[m][j]=0; }Crank_Nicolson(r,U);//調(diào)用函數(shù) printf("精確解為:u=exp(-pi^2*t)*sin(pi*x),以下打印x=0.5處一組數(shù)值解\n\n"); printf("t[n]精確解Crank_Nicolson誤差\n");//打印x=0.5處的一列點(diǎn) //for(i=1,j=1;i<=800;i=int(pow(2,j)),j++) for(i=1,j=1;i<=800;i=40*j,j++) printf("t[%3d]=%-10.8lf%-15.12lf%-15.12lf%e\n",i,t[i],u(t[i],x[10]),U[10][i],fabs(u(t[i],x[10])-U[10][i])); }voidfenjie(doublea,doubleb,doublec,doublee[m-1],doubled[m-1])//三角分解{ inti;d[0]=a; for(i=2;i<=m-1;i++) { e[i-1]=b/d[i-2]; d[i-1]=a-c*e[i-1]; }}voidhuidai(doubleg[m-1],doubley[m-1],doublee[m-1],doubled[m-1],doublec,doublex[m-1])//回代求解{ inti; for(i=2,y[0]=g[0];i<=m-1;i++) y[i-1]=g[i-1]-e[i-1]*y[i-2];//追的過程 for(i=m-2,x[m-2]=y[m-2]/d[m-2];i>=1;i--) { x[i-1]=(y[i-1]-x[i]*c)/d[i-1];//趕的過程 }}voidCrank_Nicolson(doubler,doubleU[m+1][n+1])//追趕法解決Crank_Nicolson差分問題{doublea,b,c,g[m-1],e[m-1],d[m-1],y[m-1];doublexx[m-1];inti,j;a=1+r;b=-1.0/2*r;c=b; for(j=0;j<=n;j++) {for(i=0;i<m-1;i++) { if(i==0) g[0]=(1-r)*U[1][j]-b*U[2][j]-b*(U[0][j]+U[0][j+1]); if(i==m-2) g[m-2]=(1-r)*U[m-1][j]-b*U[m-2][j]-b*(U[m][j]+U[m][j+1]); else g[i]=(1-r)*U[i+1][j]-b*(U[i+2][j]+U[i][j]); } fenjie(a,b,c,e,d);huidai(g,y,e,d,c,xx); for(i=0;i<m-1;i++) U[i+1][j+1]=xx[i]; } }doubleu(doublet,doublex)//微分方程的精確解{ doubleuu; uu=exp(pow(pi,2)*(-t))*sin(pi*x); returnuu;}附錄2：實(shí)驗(yàn)報(bào)告填寫說明1．實(shí)驗(yàn)工程名稱：要求與實(shí)驗(yàn)教學(xué)大綱一致。2．實(shí)驗(yàn)?zāi)康模耗康囊鞔_，要抓住重點(diǎn)，符合實(shí)驗(yàn)教學(xué)大綱要求。3．實(shí)驗(yàn)原理：簡要說明本實(shí)驗(yàn)工程所涉及的理論知識(shí)。4．實(shí)驗(yàn)環(huán)境：實(shí)驗(yàn)用的軟、硬件環(huán)境。5．實(shí)驗(yàn)方案〔思路、步驟和方法等〕：這是實(shí)驗(yàn)報(bào)告極其重要的內(nèi)容。概括整個(gè)實(shí)驗(yàn)過程。對(duì)于驗(yàn)證性實(shí)驗(yàn)，要寫明依據(jù)何種原理、操作方法進(jìn)行實(shí)驗(yàn)，要寫明需要經(jīng)過哪幾個(gè)步驟來實(shí)現(xiàn)