這里寫(xiě)目錄標(biāo)題

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

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

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

• 選取網(wǎng)絡(luò)；
  由于matlab數(shù)組下標(biāo)從1開(kāi)始，而離散變量的初值都是從0開(kāi)始，因此除卻初始條件，我們所求的矩陣區(qū)域?qū)?yīng)的數(shù)組下標(biāo)是 2~ N+1
  這點(diǎn)需要注意，像下文中的代碼，實(shí)際意義上是將時(shí)域劃分為了length(t)-1個(gè)區(qū)域，因?yàn)橄聵?biāo)1對(duì)應(yīng)的是初值，但只求了 2~ N
  我寫(xiě)代碼時(shí)習(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）時(shí)的循環(huán)過(guò)程則為 for i=1:m-1-1 ...

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

• 求解差分格式；

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

常微分方程用歐拉法

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

1，N+1

這里是說(shuō)算法的尺比，每個(gè)算法尺比有要求，就是空間步長(zhǎng)與時(shí)間步長(zhǎng)的比

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

求解準(zhǔn)備（對(duì)矩陣的構(gòu)造）

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



表達(dá)的式子是

matlab代碼

改成簡(jiǎn)單的邊界條件

m1=1+0.0*sin(t); %t的函數(shù)，兩個(gè)邊界條件 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的個(gè)數(shù)一樣 u(:,1)=sin(pi*x);%x的函數(shù)，初始函數(shù)，第一列 m1=1+0.0*sin(t); %t的函數(shù)，兩個(gè)邊界條件 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 主對(duì)角線，單位矩陣 %diag(,1)對(duì)角矩陣往上挪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

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

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

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





該塊內(nèi)容來(lá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求解思路

古典解、廣義解
對(duì)于PDE(偏微分方程)來(lái)說(shuō)，如果存在一個(gè)函數(shù)u uu具有所需要的各階連續(xù)偏導(dǎo)數(shù)，將它們帶入方程時(shí)能使方程成為恒等式，則稱(chēng)這個(gè)函數(shù)為該方程的解 (這種解又稱(chēng)為古典解)。

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

求解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 個(gè)分量的近似值 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時(shí)，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é)

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

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