|
|
EDA365欢迎您登录！
您需要 登录 才可以下载或查看，没有帐号？注册
x
9 |9 b5 j/ a! W% k/ V& Z& K! _MATLAB源程序代码分享：MATLAB实现偏微分方程（扩散方程）的有限差分求解法 S, T7 r& e$ F9 `' T
) E' _1 u* u/ ~; R4 a G
# M6 B1 v2 D# L, i, m+ A
* [: L# ?1 }# o" D1 P
/ `1 J/ p5 w8 ^2 h6 `%% 定义 (x,t) 平面上的网格点坐标% {8 Q& r+ `! E: ~% W- U' z: m/ q
clear;clc;close all# N0 x# u0 }3 K3 a& \$ }5 s0 N5 z
dx=0.05; % x 方向的步长+ m' Y5 J# H$ h3 A& S
dt=0.001; % t 方向的步长
9 p& i- G0 v7 B( u' G/ gx=0:dx:1; % 得到 x 的序列 (离散点 x 坐标)
+ N) D s; c% zt=0:dt:0.3; % 得到 t 的序列 (离散点 t 坐标)
8 ~7 G) K! v$ c! F$ }* Fr=dt/(dx^2); % 计算 r 的值0 @& S! K, E' A Z" s" A
/ V! A w% A" X& ]0 ^%% 设置偏微分方程的初始条件, 边界条件
/ V* z- a2 X% s! Y0 t9 c5 d" IM=length(x)-1; % 计算 x 方向的分段数 x/ z: w* c# w" S3 ~8 x
N=length(t)-1; % 计算 t 方向的分段数/ P( l( }3 w, G( {1 C9 `5 r
2 ~2 d7 m8 A4 O- V+ M7 r, r
Phi=ones(N+1,M+1); % 生成一个初始值全部为 1 的矩阵
! n- f/ S, v8 K2 D+ C) y+ IPhi(1,: )=50; % 设置初值条件: Phi(x,0)=50
* J2 z, g9 J# z& I: `Phi(2: N+1,1)=0; % 设置边界条件: Phi(0,t)=01 u6 s1 g/ E) v! _
Phi(2: N+1,M+1)=0; % 设置边界条件: Phi(1,t)=0' L+ Y) W9 ]4 B3 ]! W2 K! i
1 Z6 }- e6 X1 C
%% 根据推导出的差分方程, 计算偏微分方程的数值解6 r$ @: r6 E& m+ ?% m* m" _% d$ [+ K
for k=1:N
$ @+ K8 H: w; T2 p Y7 u. l; ~" ]/ J for i=2:M
: Q. V& w9 e" |5 }) _4 Y, n Phi(k+1,i)=(1-2*r)*Phi(k,i)+r*(Phi(k,i-1)+Phi(k,i+1));
6 I1 C7 R# H4 h" O6 |5 r7 P end
4 \& [" X) C1 t7 O3 t# Wend
2 `0 E2 H/ `# ~9 U& k0 Q, o. y% f/ a9 l: Y
%% 绘制偏微分方程的计算结果: (x, t, Phi) 的三维网格图
% r+ _; y" v6 d8 I, S! v7 Tfigure! W& k+ L/ ?3 l
set(gcf,'units','normalized','position',[0.2 0.2 0.6 0.6]); % 设置 figure 窗口的位置和尺寸
7 ~3 s7 r6 r, U$ a[x,t]=meshgrid(x,t); % 得到所有计算点的 x 坐标和 t 坐标/ |# S. C3 I3 T0 R! V. b+ ~
mesh(x,t,Phi) % 绘制 (x,t,Phi) 的三维网格图
; v; F% p* d# V- |/ W/ Pxlabel('x') 4 |( `" }4 l# @& O; ]
ylabel('t')! ^$ g5 ~$ E5 G4 r4 l) E
zlabel('\Phi(x,t)')0 m4 a R/ A; U# l- D
title('扩散方程的数值模拟'); L9 I' A% X% `3 y
view(75,50)
3 F0 Y% `$ I4 U3 r3 @; j8 Z* q7 x* K: Z4 G! b- {
|
|
/1 
发表于 2020-3-11 10:20
收藏
淘帖
支持！
反对！