這里寫目錄標題

• 微分方程數(shù)值求解——有限差分法
• matlab代碼
• • 差分法的運用（依舊是連續(xù)變量——>離散網(wǎng)格點）
• PDE求解思路
• demo1
• demo2

微分方程數(shù)值求解——有限差分法

差分方法又稱為有限差分方法或網(wǎng)格法，是求偏微分方程定解問題的數(shù)值解中應(yīng)用 最廣泛的方法之一。它的基本思想是：先對求解區(qū)域作網(wǎng)格剖分，將自變量的連續(xù)變化 區(qū)域用有限離散點（網(wǎng)格點）集代替；將問題中出現(xiàn)的連續(xù)變量的函數(shù)用定義在網(wǎng)格點 上離散變量的函數(shù)代替；通過用網(wǎng)格點上函數(shù)的差商代替導(dǎo)數(shù)，將含連續(xù)變量的偏微分 方程定解問題化成只含有限個未知數(shù)的代數(shù)方程組（稱為差分格式）。如果差分格式有 解，且當(dāng)網(wǎng)格無限變小時其解收斂于原微分方程定解問題的解，則差分格式的解就作為 原問題的近似解（數(shù)值解）。因此，用差分方法求偏微分方程定解問題一般需要解決一下問題：

• 選取網(wǎng)絡(luò)；
  由于matlab數(shù)組下標從1開始，而離散變量的初值都是從0開始，因此除卻初始條件，我們所求的矩陣區(qū)域?qū)?yīng)的數(shù)組下標是 2~ N+1
  這點需要注意，像下文中的代碼，實際意義上是將時域劃分為了length(t)-1個區(qū)域，因為下標1對應(yīng)的是初值，但只求了 2~ N
  我寫代碼時習(xí)慣于

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）時的循環(huán)過程則為 for i=1:m-1-1 ...

• 對微分方程及定解條件選擇差分近似，列出差分格式；

• 求解差分格式；

• 討論差分格式解對于微分方程解的收斂性及誤差估計

常微分方程用歐拉法

大部分，下標位置，上標時間
有熱源或者是擴散源，連續(xù)的、分成很多份、離散地取值，有限多個數(shù)來近似一個函數(shù)

1，N+1

這里是說算法的尺比，每個算法尺比有要求，就是空間步長與時間步長的比

初始算出 1，1~ N+1，可得到2，2~N，通過 邊界條件，可得到 2，1 ~ N+1,繼續(xù)向上推

求解準備（對矩陣的構(gòu)造）

左側(cè)矩陣是u關(guān)于x的二階導(dǎo)，但只是先后對于 2~ N的二階導(dǎo)，要求得 1 ~ N+1 的二階導(dǎo)，還需要借助邊界條件



表達的式子是

matlab代碼

改成簡單的邊界條件

m1=1+0.0*sin(t); %t的函數(shù)，兩個邊界條件 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)); %行數(shù)和x的個數(shù)一樣 u(:,1)=sin(pi*x);%x的函數(shù)，初始函數(shù)，第一列 m1=1+0.0*sin(t); %t的函數(shù)，兩個邊界條件 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)-1u(:,n+1)=u(:,n)+a^2*dt/dx^2*A*u(:,n);% a^2*(u關(guān)于x的二階導(dǎo))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

差分法的運用（依舊是連續(xù)變量——>離散網(wǎng)格點）

由于向前差分有誤差，如果我們進行兩次向前差分的話，計算的誤差可能會增大，因此，第二次偏導(dǎo)我們選擇向后差分。即我們混合向前差分、向后差分來近似代替兩次偏導(dǎo)。

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





該塊內(nèi)容來自博文

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;% 總列數(shù) n=4;% 總行數(shù) 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(偏微分方程)來說，如果存在一個函數(shù)u uu具有所需要的各階連續(xù)偏導(dǎo)數(shù)，將它們帶入方程時能使方程成為恒等式，則稱這個函數(shù)為該方程的解 (這種解又稱為古典解)。

用一個充分光滑的初值函數(shù)序列來逼近不夠光滑的初值函數(shù)，前者所對應(yīng)的解序列的極限就定義為后者所確定的解，稱為問題的廣義解。

求解ODE思路
求解常微分方程的辦法，先求出方程的通解，再用定解條件去確定任意常數(shù)。現(xiàn)在，如能找出主方程的通解，再利用定解條件去確定任意函數(shù)。

求解PDE思路
求出PDE滿足邊界條件的足夠數(shù)目的特解，再利用疊加原理，使之滿足初始條件，從而得到混合方程的解。

工具箱求解

demo1

參考博文

demo2

%主函數(shù) function main clc,clear m=0; x=linspace(0,1,20); % 方程區(qū)間為(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文檔

總結(jié)

以上是生活随笔為你收集整理的差分、偏微分方程的解法的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。