差分、偏微分方程的解法

微分方程数值求解——有限差分法

差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用 最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化 区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点 上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分 方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有 解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为 原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决一下问题：

• 选取网络；
  由于matlab数组下标从1开始，而离散变量的初值都是从0开始，因此除却初始条件，我们所求的矩阵区域对应的数组下标是 2~ N+1
  这点需要注意，像下文中的代码，实际意义上是将时域划分为了length(t)-1个区域，因为下标1对应的是初值，但只求了 2~ N
  我写代码时习惯于

dx=0.02;
x=0:dx:1;
dt=0.0001;%0.0001
t=0:dt:1;
n=length(x)+1,m=length(t)+1
u=zeros(n,m);
%在求解U（x，t）时的循环过程则为
for i=1:m-1-1
...

• 对微分方程及定解条件选择差分近似，列出差分格式；

• 求解差分格式；

• 讨论差分格式解对于微分方程解的收敛性及误差估计

常微分方程用欧拉法

大部分，下标位置，上标时间
有热源或者是扩散源，连续的、分成很多份、离散地取值，有限多个数来近似一个函数

1，N+1

这里是说算法的尺比，每个算法尺比有要求，就是空间步长与时间步长的比

初始算出 1，1~ N+1，可得到2，2~N，通过 边界条件，可得到 2，1 ~ N+1,继续向上推

求解准备（对矩阵的构造）

左侧矩阵是u关于x的二阶导，但只是先后对于 2~ N的二阶导，要求得 1 ~ N+1 的二阶导，还需要借助边界条件



表达的式子是

matlab代码

改成简单的边界条件

m1=1+0.0*sin(t); %t的函数，两个边界条件
m2=2-0.0*sin(10*t);

clc,clear
a=1;
dx=0.02;
x=0:dx:1;
dt=0.0001;%0.0001
t=0:dt:1;
u=zeros(length(x),length(t));
%行数和x的个数一样
u(:,1)=sin(pi*x);%x的函数，初始函数，第一列
m1=1+0.0*sin(t); %t的函数，两个边界条件
m2=2-0.0*sin(10*t);
A=-2*eye(length(x))+diag(ones(1,length(x)-1),1)++diag(ones(1,length(x)-1),-1);
%eye 主对角线，单位矩阵
%diag(,1)对角矩阵往上挪1
for n=1:length(t)-1
    u(:,n+1)=u(:,n)+a^2*dt/dx^2*A*u(:,n);% a^2*(u关于x的二阶导)
    u(1,n+1)=m1(n+1); %边界条件
    u(end,n+1)=m2(n+1);
    
end

%plot(x,u(:,end),'-bp')

  [X,T]=meshgrid(t,x);
    surf(X,T,u)
    shading interp  

差分法的运用（依旧是连续变量——>离散网格点）

由于向前差分有误差，如果我们进行两次向前差分的话，计算的误差可能会增大，因此，第二次偏导我们选择向后差分。即我们混合向前差分、向后差分来近似代替两次偏导。

因此，第二次我们用向后差分





该块内容来自博文

clc;
clear;
f1 = @(x)2 * log(1+x);
f2 = @(x)log((1+x).^2+1);
f3 = @(y)log(1+y.^2);
f4=  @(y)log(4+y.^2);
u=zeros(4);
m=4;% 总列数
n=4;% 总行数
h=1/3;
u(1,1:m)=feval(f3,0:h:(m-1)*h)';
u(n,1:m)=feval(f4,0:h:(m-1)*h)';
u(1:n,1)=feval(f1,0:h:(n-1)*h);
u(1:n,m)=feval(f2,0:h:(n-1)*h);
b = -[u(2,1)+u(1,2);u(4,2)+u(3,1);u(2,4)+u(1,3);u(3,4)+u(4,3)];
a = [-4,1,1,0;1,-4,0,1;1,0,-4,1;0,1,1,-4];
x=a\b;

PDE求解思路

古典解、广义解
对于PDE(偏微分方程)来说，如果存在一个函数u uu具有所需要的各阶连续偏导数，将它们带入方程时能使方程成为恒等式，则称这个函数为该方程的解 (这种解又称为古典解)。

用一个充分光滑的初值函数序列来逼近不够光滑的初值函数，前者所对应的解序列的极限就定义为后者所确定的解，称为问题的广义解。

求解ODE思路
求解常微分方程的办法，先求出方程的通解，再用定解条件去确定任意常数。现在，如能找出主方程的通解，再利用定解条件去确定任意函数。

求解PDE思路
求出PDE满足边界条件的足够数目的特解，再利用叠加原理，使之满足初始条件，从而得到混合方程的解。

工具箱求解

demo1

参考博文

demo2

%主函数
function main
clc,clear
m=0;
x=linspace(0,1,20); % 方程区间为(0,1)
t=linspace(0,2,10); % t 的范围可以随取，只需要大于0即可
sol=pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % 10x20的矩阵与t*x的维度一致
u = sol(:,:,1); %解向量 u 的第 1 个分量的近似值
surf(x,t,u) 
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

figure
plot(x,u(end,:))% t=2时，u随x的变化曲线
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,2)')

end
%PDE方程 
function [c,f,s]=pdex1pde(x,t,u,DuDx)
    c=pi^2;
    f=DuDx;
    s=0;
 end
%初始条件格式
function uo=pdex1ic(x)
    uo=sin(pi*x);
end
%边界条件
function [pl,ql,pr,qr]=pdex1bc(x1,u1,xr,ur,t)
    pl=u1; 
    ql=0; 
    pr=pi*exp(-t); 
    qr=1; 
    end

参考博文
matlab文档