1 問題描述

Chenglin Li：數值計算（三）matlab求解一般的偏微分方程組?zhuanlan.zhihu.com

• 因為給出的邊界條件包含導數，因此需要同時考慮前向差分和后向差分；
• 遍歷循環，先計算每個坐標的時間節點，或者先計算每個時間節點的坐標，結果不一樣；

2 計算程序（太復雜未完成，給出思路）

function [pa, U]=pde201105() %{ 程序功能： 1、熱傳導PDE方程組的差分求解 2、0<x<a, 0<t<b. 3、date:2020.11.20 %}clear,clc, close all %參數聲明pa.a=1 ; %0<x<apa.b=0.2;%1/3;%0.2 ; %0<t<bpa.h=0.1;pa.k=0.02;%1/30; %0.02;pa.n=pa.a/pa.h+1 ; %x坐標節點數pa.m=pa.b/pa.k+1; %t時間節點數pa.c1=0.024 ;pa.c2=0.170 ;pa.r1=pa.c1/pa.k ;pa.r2=pa.c2/pa.k ;pa.s1=2-1/pa.r1 ;pa.s2=2-1/pa.r2 ;u=zeros(pa.n, pa.m, 2) ; %x-t%-------------------------%初值約束for i=2:pa.n-1y=fx( pa.h* (i-1) );u(i, 1, 1) =y(1);u(i,1, 2) =y(2) ;end%邊值約束1-常數ufor j=2: pa.my=ga( pa.k*(j-1) ) ;u(1,j,1)=y(1) ;y=gb( pa.k*(j-1) ) ;u(pa.n, j, 2) =y(2) ;end%邊值約束2-導數dufor j=2: pa.my=pa.h*dga( pa.k*(j-1) );u(2, j, 1)=y(1)+u(1,j,1);y=pa.h*dgb( pa.k*(j-1) );u(pa.n-1, j, 2)=y(2)-u(pa.n, j, 2) ;end%依次計算每一時刻的節點值%正序差分for i=3: pa.n-1 %先計算行--位移for j=2: pa.m %次計算列 --時間u(i, j,1)=pa.s1*u(i-1, j, 1)+1/pa.r1*( u(i-1, j+1,1)+ pa.k*Fx( u(i-1, j,1)-u(i-1, j, 2 ) ) )-u(i-2, j, 1) ; %正序插值endend%倒序差分for i=pa.n-2: -1 :2for j=2: pa.mu(i, j,2)=pa.s2*u(i-1, j, 2)+1/pa.r2*( u(i-1, j+1,2)+ pa.k*Fx( u(i-1, j,1)-u(i-1, j, 2) ) ) -u(i-2, j, 1) ; %倒序插值endend%繪制圖形pa.x=0: pa.h: pa.a;pa.t=0: pa.k: pa.b ;U(:,:,1)=u(:,:,1)' ;U(:,:,2)=u(:,:,2)';[X,T]=meshgrid(pa.x, pa.t) ;figure(1)mesh(X,T,U(:,:,1))xlabel('x')ylabel('t')zlabel('u1')figure(2)mesh(X,T,U(:,:,2))xlabel('x')ylabel('t')zlabel('u1')end%右端項 function y=Fx( x )y=exp(5.73*x) -exp(-11.46*x) ;end %初始條件：u(x,0)=[f1(x),f2(x) ] function y=fx(x) % y=sin(pi*x)+sin(2*pi*x) ; % y=4*x-4*x^2; % y=sin(pi*x)+sin(3*pi*x);y(1)=1 ; y(2)=0 ;end%u(0,t)=g1(t) function y=ga(t)y(1)=0;y(2)=NaN;end function y=dga(t)y(1)=0;y(2)=NaN;end %u(a,t)=g2(t) function y=gb(t)y(1)=NaN;y(2)=1;endfunction y=dgb(t)y(1)=NaN ;y(2)=0 ; end

——2020.11.20——

總結

以上是生活随笔為你收集整理的计算差分方程的收敛点_数值计算（五十九）热传导方程组的差分数值求解的全部內容，希望文章能夠幫你解決所遇到的問題。

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