日韩av黄I国产麻豆传媒I国产91av视频在线观看I日韩一区二区三区在线看I美女国产在线I麻豆视频国产在线观看I成人黄色短片

歡迎訪問 生活随笔!

生活随笔

當前位置: 首頁 >

差分、偏微分方程的解法

發布時間:2024/5/8 94 豆豆
生活随笔 收集整理的這篇文章主要介紹了 差分、偏微分方程的解法 小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

這里寫目錄標題

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

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

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

  • 選取網絡;
    由于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)-1u(:,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文檔

總結

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

如果覺得生活随笔網站內容還不錯,歡迎將生活随笔推薦給好友。