偏微分方程 数值解

偏微分方程数值解法学生姓名：学号：学院：专业年级：设计题目：一维热传导方程几种差分格式2012年07月一维热传导方程定解条件：1．初值问题（Cauchy问题）2．初边值问题（混合问题）物理意义：有限长细杆的温度（两端温度不一定为0，一般为，这里取成0为方便计算）方程中的条件为第一边值条件，也可有第二、三类边值条件：这里的系数应满足一定的条件（见书）。它们实际上包含了第一、第二、第三类边界条件。假定，在区域上充分光滑，则问题的解存在且唯一。现考虑边值问题的差分格式。对区域进行分割：用平行直线族分割，为节点，记为，称为空间步长，称为时间步长。记用表示定义在网点的网函数，，。用不同的差商代替方程中的偏微商，可得差分格式。一.向前差分格式（古典显格式）两边同乘以，记（称为网比）则上式可写成便于计算的形式：（层在左边，层在右边）取，依次可算出各层上的值。此格式不必求解方程组，直接可算，称显格式。，可得截断误差阶为二向后差分格式（古典隐格式）同样可化成：（网比）每一方程有3个未知量，要解三对角阵的线性方程组，计算量比显格式大，系数矩阵严格对角占优，可证明其非奇异，此方程组有唯一解。截断误差。三六点对称格式（Crank-Nicolson格式）将向前差分格式和向后差分格式作算术平均，得：边值条件的处理同上。记网比，化成：每一层也要计算一个线性代数方程组的解，也是隐格式，系数也是对角占优，有唯一解。截断误差阶：。注：将向前差分格式乘以，向后差分格式格式乘以相加，所得格式称加权六点格式。当时，即为六点对称格式。对加权六点格式（设）其中，书上P.94-95用泰勒展开的方法推导了它的截断误差：则当时，截断误差的阶最高（）。四Richardson格式：用中心差商代替对的微商，得：记，化为。这是一个三层的差分格式。已知，还需用其他两层格式求出，才能用此格式求出。它是显格式，截断误差的阶是。还可以造出一些其他的格式。如在上面的Richardson格式中，令代入，可得DuFortFrankel格式。哪种格式更好？是否有实用价值？为此要考虑：（1）计算简单：显格式最简单，隐格式要解方程组，可用追赶法。但比显格式麻烦。（2）收敛性和收敛速度：当网比固定时，步长（隐含），差分解应收敛到真解，即，其中。收敛速度与截断误差阶有关。应比在同样的网比时精度更高一些。（3）稳定性：计算中产生的误差在以后的计算中能否被控制？不能控制是不稳定的，不能应用。定义：，（与，无关），当时，有成立。其中是由初始向量产生的新的解，是差分方程的真解。则称差分格式是稳定的。对Richardson格式，考察其稳定性。设为计算产生的误差，设无误差，且层前和本层中其它点的计算均无误差，而计算过程中也不再产生新的误差。取，则满足：设时，中间点的计算有误差，即，。则由上式，得，，……，等等，绝对值越来越大。故此格式不稳定。实际上对任意，此格式均不稳定。对于向前差分格式，当时，误差方程为。取，，则，，……。随增加误差越来越小。五数值例子例1令f(x)=0和a=1，可求得u（x,t）一个解析解为u(x,t)=exp(x+t)。用Richardson格式验证数值结果如下：请输入n的值(输入0结束程序):5请输入m的值(输入0结束程序):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]陆金甫，关治编.偏微分方程数值解法[M].北京：清华大学出版社，2003：35-137.[2]南京大学数学系.计算数学专业编，偏微分方程数值解法[M].北京：科学出版社，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;//将空间区间【0，1】分为n等份；时间区间【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()//一维热传导方程的Richardson格式{ intk,i; doubleh,t,r; doublepi=3.1415627; printf("请输入n的值(输入0结束程序):\n"); if(scanf("%d",&n))printf("请输入m的值(输入0结束程序):\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]