偏微分方程 數(shù)值解

偏微分方程數(shù)值解法學(xué)生姓名：學(xué)號：學(xué)院：專業(yè)年級：設(shè)計題目：一維熱傳導(dǎo)方程幾種差分格式2012年07月一維熱傳導(dǎo)方程定解條件：1．初值問題（Cauchy問題）2．初邊值問題（混合問題）物理意義：有限長細桿的溫度（兩端溫度不一定為0，一般為，這里取成0為方便計算）方程中的條件為第一邊值條件，也可有第二、三類邊值條件：這里的系數(shù)應(yīng)滿足一定的條件（見書）。它們實際上包含了第一、第二、第三類邊界條件。假定，在區(qū)域上充分光滑，則問題的解存在且唯一。現(xiàn)考慮邊值問題的差分格式。對區(qū)域進行分割：用平行直線族分割，為節(jié)點，記為，稱為空間步長，稱為時間步長。記用表示定義在網(wǎng)點的網(wǎng)函數(shù)，，。用不同的差商代替方程中的偏微商，可得差分格式。一.向前差分格式（古典顯格式）兩邊同乘以，記（稱為網(wǎng)比）則上式可寫成便于計算的形式：（層在左邊，層在右邊）取，依次可算出各層上的值。此格式不必求解方程組，直接可算，稱顯格式。，可得截斷誤差階為二向后差分格式（古典隱格式）同樣可化成：（網(wǎng)比）每一方程有3個未知量，要解三對角陣的線性方程組，計算量比顯格式大，系數(shù)矩陣嚴(yán)格對角占優(yōu)，可證明其非奇異，此方程組有唯一解。截斷誤差。三六點對稱格式（Crank-Nicolson格式）將向前差分格式和向后差分格式作算術(shù)平均，得：邊值條件的處理同上。記網(wǎng)比，化成：每一層也要計算一個線性代數(shù)方程組的解，也是隱格式，系數(shù)也是對角占優(yōu)，有唯一解。截斷誤差階：。注：將向前差分格式乘以，向后差分格式格式乘以相加，所得格式稱加權(quán)六點格式。當(dāng)時，即為六點對稱格式。對加權(quán)六點格式（設(shè)）其中，書上P.94-95用泰勒展開的方法推導(dǎo)了它的截斷誤差：則當(dāng)時，截斷誤差的階最高（）。四Richardson格式：用中心差商代替對的微商，得：記，化為。這是一個三層的差分格式。已知，還需用其他兩層格式求出，才能用此格式求出。它是顯格式，截斷誤差的階是。還可以造出一些其他的格式。如在上面的Richardson格式中，令代入，可得DuFortFrankel格式。哪種格式更好？是否有實用價值？為此要考慮：（1）計算簡單：顯格式最簡單，隱格式要解方程組，可用追趕法。但比顯格式麻煩。（2）收斂性和收斂速度：當(dāng)網(wǎng)比固定時，步長（隱含），差分解應(yīng)收斂到真解，即，其中。收斂速度與截斷誤差階有關(guān)。應(yīng)比在同樣的網(wǎng)比時精度更高一些。（3）穩(wěn)定性：計算中產(chǎn)生的誤差在以后的計算中能否被控制？不能控制是不穩(wěn)定的，不能應(yīng)用。定義：，（與，無關(guān)），當(dāng)時，有成立。其中是由初始向量產(chǎn)生的新的解，是差分方程的真解。則稱差分格式是穩(wěn)定的。對Richardson格式，考察其穩(wěn)定性。設(shè)為計算產(chǎn)生的誤差，設(shè)無誤差，且層前和本層中其它點的計算均無誤差，而計算過程中也不再產(chǎn)生新的誤差。取，則滿足：設(shè)時，中間點的計算有誤差，即，。則由上式，得，，……，等等，絕對值越來越大。故此格式不穩(wěn)定。實際上對任意，此格式均不穩(wěn)定。對于向前差分格式，當(dāng)時，誤差方程為。取，，則，，……。隨增加誤差越來越小。五數(shù)值例子例1令f(x)=0和a=1，可求得u（x,t）一個解析解為u(x,t)=exp(x+t)。用Richardson格式驗證數(shù)值結(jié)果如下：請輸入n的值(輸入0結(jié)束程序):5請輸入m的值(輸入0結(jié)束程序):5xjtk真實值x[i][k]近似值u[i][k]誤差err[i][k]0.1666670.1666671.3956121.3960800.0004680.3333330.1666671.6487211.6494810.0007600.5000000.1666671.9477341.9486490.0009150.6666670.1666672.3009762.3018880.0009120.8333330.1666672.7182822.7189490.0006670.1666670.3333331.6487211.6455400.0031820.3333330.3333331.9477341.9448060.0029280.5000000.3333332.3009762.2975830.0033930.6666670.3333332.7182822.7136000.0046820.8333330.3333333.2112713.2040990.0071710.1666670.5000001.9477341.9881450.0404110.3333330.5000002.3009762.2916210.0093550.5000000.5000002.7182822.7075140.0107680.6666670.5000003.2112713.1956850.0155860.8333330.5000003.7936683.9077790.1141110.1666670.6666672.3009761.2141561.0868200.3333330.6666672.7182823.2938220.5755410.5000000.6666673.2112713.1649070.0463640.6666670.6666673.7936685.4006921.6070240.8333330.6666674.4816891.5458782.9358110.1666670.8333332.71828235.74707433.0287920.3333330.8333333.211271-24.21136127.4226310.5000000.8333333.79366831.08392727.2902590.6666670.8333334.481689-69.89150974.3731980.8333330.8333335.29449095.14889189.854401六．參考文獻[1]陸金甫，關(guān)治編.偏微分方程數(shù)值解法[M].北京：清華大學(xué)出版社，2003：35-137.[2]南京大學(xué)數(shù)學(xué)系.計算數(shù)學(xué)專業(yè)編，偏微分方程數(shù)值解法[M].北京：科學(xué)出版社，1979：10-11.附錄程序源代碼：#include<stdio.h>#include<stdlib.h>#include<math.h>#defineMax_N1000double u[Max_N][Max_N],b[Max_N],a[Max_N],c[Max_N],f[Max_N],err[Max_N][Max_N],x[Max_N][Max_N],y[Max_N],beta[Max_N],Err[Max_N];intn,m;//將空間區(qū)間【0，1】分為n等份；時間區(qū)間【0，1】分為m等份voidcatchup(){ inti; beta[1]=c[1]/b[1]; for(i=2;i<n;i++) beta[i]=c[i]/(b[i]-a[i]*beta[i-1]); y[1]=f[1]/b[1]; for(i=2;i<=n;i++) y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]); u[n][1]=y[n]; for(i=n-1;i>0;i--) u[i][1]=y[i]-beta[i]*u[i+1][1];}intmain()//一維熱傳導(dǎo)方程的Richardson格式{ intk,i; doubleh,t,r; doublepi=3.1415627; printf("請輸入n的值(輸入0結(jié)束程序):\n"); if(scanf("%d",&n))printf("請輸入m的值(輸入0結(jié)束程序):\n"); while(scanf("%d",&m)&&m&&n)//u(x,t)=exp(x+t),u(x,0)=exp(x),f(x)=0,x屬于[0,1],t屬于[0,1],a=1. { h=1.0/(n+1); t=1.0/(m+1); r=t/(h*h); for(i=0;i<=n+1;i++)//初值條件 { u[i][0]=exp(i*h); } for(k=0;k<=m+1;k++)//邊值條件 { u[0][k]=exp(k*t); u[n+1][k]=exp((n+1)*h+k*t); } printf("xjtk真實值x[i][k]近似值u[i][k]誤差err[i][k]\n"); b[1]=1+r; c[1]=-r/2; a[n]=-r/2; b[n]=1+r; f[1]=r/2*u[2][0]+(1-r)*u[1][0]+r/2+r/2*u[0][1]; f[n]=r/2*u[n+1][0]+(1-r)*u[n][0]+r/2*u[n-1][0]+r/2*u[n+1][1]; for(i=2;i<n;i++) { b[i]=1+r; a[i]=-r/2; c[i]=-r/2; f[i]=r/2*u[i+1][0]+(1-r)*u[i][0]+r/2*u[i-1][0]; } catchup(); for(k=2;k<=m;k++) { for(intj=1;j<=n;j++) { u[j][k]=2*r*(u[j+1][k-1]-2*u[j][k-1]+u[j-1][k-1])+u[j][k-2]; } } for(k=1;k<=m;k++) { for(i=1;i<=n;i++) { x[i][k]=exp(i*h+k*t); err[i]