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——
