问题内容

用显格式，向隐格式Crank-Nicolson格式来求解，取h=0.1，r=为0.1进行计算，后并分析误差。（τ\tauτ为t步长，h为x步长，为m位置，k时刻的温度）

算法求解

显格式（explicit scheme）求解

参考课本可得显格式（explicit scheme）的迭代过程：
(u1k+1u2k+1⋮um−2k+1um−1k+1)=(1−2rrr1−2rr⋱⋱⋱r1−2rrr1−2r)(u1k+ru0ku2k⋮um−2kum−1k+rumk)\left(\begin{array}{c}u_{1}^{k+1} \\ u_{2}^{k+1} \\ \vdots \\ u_{m-2}^{k+1} \\ u_{m-1}^{k+1}\end{array}\right)=\left(\begin{array}{ccccc}1-2 r & r & & & \\ r & 1-2 r & r & & \\ & \ddots & \ddots & \ddots & \\ & & & r & 1-2 r & r \\ & & & & r & 1-2 r\end{array}\right)\left(\begin{array}{c}u_{1}^{k}+ru_{0}^k \\ u_{2}^{k} \\ \vdots \\ u_{m-2}^{k} \\ u_{m-1}^{k}+ru_{m}^k\end{array}\right)⎝⎜⎜⎜⎜⎜⎛u1k+1u2k+1⋮um−2k+1um−1k+1⎠⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎛1−2rrr1−2r⋱r⋱⋱r1−2rrr1−2r⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1k+ru0ku2k⋮um−2kum−1k+rumk⎠⎟⎟⎟⎟⎟⎞
实现程序：

%参考课本
%slove program by explicit scheme
%u(j,n+1)=u(j,n)+v*(u(j+1,n)-2*u(j,n)+u(j-1,n))
xl=0;
xr=1;
j=10;
dx=(xr-xl)/j;%步长
tf=0.1;
Nt=50;
dt=tf/Nt;
mu=dt/(dx)^2;
%make sure dt satisy stability condition
if mu>0.5error('mu shuold<0.5!')
end
%initial condition
x=xl:dx:xr;%grind point
f=sin(pi*x)+sin(2*pi*x);
%store the solution at all grid points for all time steps
u=zeros(j+1,Nt);
u_ture=zeros(j+1,Nt);
%find the approximate solution at each time step
for n=1:Ntt=n*dt;%boundary condition at left sidegl=exp(-1*pi*pi*t).*sin(pi*xl)+exp(-4*pi*pi*t).*sin(2*pi*xl);%boundary condition at right sidegr=exp(-1*pi*pi*t).*sin(pi*xr)+exp(-4*pi*pi*t).*sin(2*pi*xr);if n==1for i=2:ju(i,n)=f(i)+mu*(f(j+1)-2*f(j)+f(j-1));endu(1,n)=gl;u(j+1)=gr;elsefor i=2:ju(i,n)=u(i,n-1)+mu*(u(i+1,n-1)-2*u(i,n-1)+u(i-1,n-1));endu(1,n)=gl;u(j+1)=gr;end%calculate the analytic solutionfor i=1:j+1xi=xl+(i-1)*dx;u_ture(i,n)=exp(-1*pi*pi*t).*sin(pi*xi)+exp(-4*pi*pi*t).*sin(2*pi*xi);endend
%plot the result
tt=dt:dt:Nt*dt;
figure(1)
colormap(gray);
surf(x,tt,u');
xlabel('x');
ylabel('t');
zlabel('u');
title('explicit scheme')%polt the analytic result
figure(2)
colormap(jet);
surf(x,tt,u_ture');
xlabel('x');
ylabel('t');
zlabel('u');
title('analytic solution')

隐格式（implicit scheme）求解

参考课本可得隐格式（implicit scheme）的迭代过程：
(1+2r−r−r1+2r−r⋱⋱⋱−r1+2r−r−r1+2r)(u1ku2k⋮um−2kum−1k)\left(\begin{array}{ccccc}1+2 r & -r & & & \\ -r & 1+2 r & -r & & \\ & \ddots & \ddots & \ddots & \\ & & -r & 1+2 r & -r \\ & & & -r & 1+2 r\end{array}\right)\left(\begin{array}{c}u_{1}^{k} \\ u_{2}^{k} \\ \vdots \\ u_{m-2}^{k} \\ u_{m-1}^{k}\end{array}\right)⎝⎜⎜⎜⎜⎛1+2r−r−r1+2r⋱−r⋱−r⋱1+2r−r−r1+2r⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1ku2k⋮um−2kum−1k⎠⎟⎟⎟⎟⎟⎞
=(u1k−1+ru0ku2k−1⋮um−2k−1um−1k−1+rumk)=\left(\begin{array}{c}u_{1}^{k-1}+r u_{0}^{k} \\ u_{2}^{k-1} \\ \vdots \\ u_{m-2}^{k-1} \\ u_{m-1}^{k-1}+r u_{m}^{k}\end{array}\right)=⎝⎜⎜⎜⎜⎜⎛u1k−1+ru0ku2k−1⋮um−2k−1um−1k−1+rumk⎠⎟⎟⎟⎟⎟⎞

参照上算法

Crank_Nicolsm格式求解

参考课本可得Crank_Nicolsm格式的迭代过程：
(1+r−r2−r21+r−r2⋱⋱⋱−r21+r−r2−r21+r)(u1k+1u2k+1⋮um−2k+1um−1k+1)\left(\begin{array}{ccccc}1+r & -\frac{r}{2} & & & \\ -\frac{r}{2} & 1+r & -\frac{r}{2} & & \\ & \ddots & \ddots & \ddots & \\ & & -\frac{r}{2} & 1+r & -\frac{r}{2} \\ & & & -\frac{r}{2} & 1+r\end{array}\right)\left(\begin{array}{c}u_{1}^{k+1} \\ u_{2}^{k+1} \\ \vdots \\ u_{m-2}^{k+1} \\ u_{m-1}^{k+1}\end{array}\right)⎝⎜⎜⎜⎜⎛1+r−2r−2r1+r⋱−2r⋱−2r⋱1+r−2r−2r1+r⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1k+1u2k+1⋮um−2k+1um−1k+1⎠⎟⎟⎟⎟⎟⎞
=(1−rr2r21−rr2⋱⋱⋱r21−rr2r21−r)(u1k+r2(u0k+u0k+1)u2k⋮um−2kum−1k+r2(umk+umk+1))=\left(\begin{array}{ccccc}1-r & \frac{r}{2} & & & \\ \frac{r}{2} & 1-r & \frac{r}{2} & & \\ & \ddots & \ddots & \ddots & \\ & & & \frac{r}{2} & 1-r & \frac{r}{2} \\ & & & \frac{r}{2} & 1-r\end{array}\right)\left(\begin{array}{c}u_{1}^{k}+\frac{r}{2}\left(u_{0}^{k}+u_{0}^{k+1}\right) \\ u_{2}^{k} \\ \vdots \\ u_{m-2}^{k} \\ u_{m-1}^{k}+\frac{r}{2}\left(u_{m}^{k}+u_{m}^{k+1}\right)\end{array}\right)=⎝⎜⎜⎜⎜⎛1−r2r2r1−r⋱2r⋱⋱2r2r1−r1−r2r⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛u1k+2r(u0k+u0k+1)u2k⋮um−2kum−1k+2r(umk+umk+1)⎠⎟⎟⎟⎟⎟⎞
程序

参照上程序

实验结果

三格式的迭代求解结果

误差分析

本次实验以CN格式分析为例，按上述计算思想计算得误差数据汇总的得到表1：（E∞=max⁡∣u(xi,tk)−uik∣\boldsymbol{E}_{\infty}=\max \left|\boldsymbol{u}\left(\boldsymbol{x}_{i}, \boldsymbol{t}_{k}\right)-\boldsymbol{u}_{i}^{k}\right|E∞=max∣∣u(xi,tk)−uik∣∣
）

实验总结

本次使用三种不同的差分格式：Crank_Nicolsm格式、显格式（explicit scheme）和隐格式（implicit scheme）对一个详细的一维抛物型方程进行求解，通过对上述格式进行误差分析，可以得到CN格式的求解误差最小，同时改变步长计算误差，得到CN格式下的误差收敛速度为：
。但是实验在三种格式的收敛速度的进一步比较上，限于时间精力，没有进一步展开，有待进一步研究提高。

一维抛物型方程的差分解法相关推荐

抛物型方程的差分解法matlab,抛物型方程的差分解法
? 0 时 2 为前差方程,当 ? ?1 时为后差方程.用控制体积法构造差分方程总是守恒型差分方程. (4) 积分方法采用积分方法构造差分方程基本思想是把微分...... 第7卷第4期2008年12 ...
椭圆形方程的差分解法
最近看<热工流体数值计算>(俞冀阳,2012)中讲了一些典型的偏微分方程求解,按照其差分计算的方法,采用Java语言和Python语言进行了式2-8的求解.数学模型如下: 差分计算时,收敛 ...
隐方程求解一维抛物型方程（热传导方程）
记录一下数值计算课大作业设时间t和空间x 离散:xi_{i}i=ih , tn_{n}n=nΔ\DeltaΔ*t (i = 0,1,-I , n = 0,1,-) 设 U(xi_{i}i,tn ...
matlab抛物偏微分方程,抛物型方程差分求解跪求MATLAB解抛物型偏微分方程的程序...
为什么抛物线方程与圆方程联立不能使用韦达定理很容易了解到,抛物线和圆的交点均在X轴上方\"其实这时你应该注意到一点就是,这两个交点的纵坐标是相等的,所以其实对应的是一个y值,也就是你列的一 ...
抛物型方程的有限差分 C语言程序,抛物型方程有限差分方法的应用 - 报告.doc
抛物型方程有限差分方法的应用 - 报告 2015 年秋季学期研究生课程考核 (读书报告.研究报告) 考核科目: 偏微分方程数值解法学生所在院(系): 理学院数学系学生所在学科: 数学学生姓 ...
PDE抛物型方程数值解法总结与例题分析
抛物型方程例题及解答例题: 构造抛物型方程 {∂u∂t=∂∂x(x∂u∂x),0.5<x<1,0<t⩽T,u(x,0)=φ(x),0.5⩽x⩽1,u(0.5,t)=0,∂u∂x( ...
微分方程近似解法归纳——差分解法
文章目录前言一.差分方程的概念二.差分格式的构造三.差分方程求解 1.基本原理: (1)先考虑一维情形 (2)二维情形第一类边界条件第二.三类边界条件(以下内容摘自教材) 总结前言数学 ...
非线性微分方程有限差分解法
一般的, 对于线性偏微分方程其有限差分解法比较容易编写. 对于非线性微分方程的数值求解而言, 一般的, 如果采用有限差分的方法进行离散, 那么离散后方程中就会出现非线性项. 此时可以将非线性项看作是新 ...
matlab 差分解微分,基于MATLAB的偏微分方程差分解法
<基于MATLAB的偏微分方程差分解法>由会员分享,可在线阅读,更多相关<基于MATLAB的偏微分方程差分解法(12页珍藏版)>请在人人文库网上搜索. 1.基于MATLAB的偏 ...

一维抛物型方程的差分解法