Matlab數值計算最簡單的一個例子——指數衰減


放射性衰變是指數衰減的典型例子。另外還有化學反應某反應物的減少,RC電路電流的減小,大氣壓隨海拔高度的減小等。

指數衰減的方程:

\begin{equation}
\frac{dN(t)}{dt}=-\frac{N(t)}{\tau}
\label{eq1}
\end{equation}

其中,\(N(t)\)\(t\)時刻的物理量\(N\),對於放射性衰變,\(N\)就是未衰變的原子核數目。\(\tau\)為時間常數。

方程\eqref{eq1}有解析解:

\[N(t)=N(0)\exp(-t/\tau) \]

這個解可以通過Matlab符號計算求得:

  dsolve('DN=-N/tau')
  ans =
  C3*exp(-t/tau)

數值求解方程\eqref{eq1},可用歐拉格式將方程離散化。

\[t_i=(i-1) \Delta t,\quad i=1,2,\dots,\mathrm{npoints} \]

\[\frac{dN(t)}{dt}\approx\frac{N(t)-N(t-\Delta t)}{\Delta t} \]

將以上兩式帶入方程\eqref{eq1},得離散之后的方程:

\[N(t_{i+1})=N(t_i)-N(t_i)\frac{\Delta t}{\tau} \]

代入初始條件,即可得解。

下面寫個Matlab 腳本文件,重復出Computational Physics_Giordano 2nd Edition的圖1.1,pp11

%
% Exponent decay
% 'Computational Physics' book by N Giordano and H Nakanishi
% Section 1.2 p2
% Solve the Equation dN/dt = -N/tau
% by Joyful Physics Blog
% ------------------------------------------------------------

N_nuclei_initial = 100; %initial number of nuclei
npoints = 101; % Discretize time into 100 intervals
dt = 0.05; % set time step 
tau=1; % set time constant
N_nuclei = zeros(npoints,1); % initializes N_nuclei, a vector of dimension npoints X 1,to being all zeros
time = zeros(npoints,1); % this initializes the vector time to being all zeros
N_nuclei(1) = N_nuclei_initial; % the initial condition, first entry in the vector N_nuclei is N_nuclei_initial
time(1) = 0; %Initialise time
for step=1:npoints-1 % loop over the timesteps and calculate the numerical solution
N_nuclei(step+1) = N_nuclei(step) - (N_nuclei(step)/tau)*dt;
time(step+1) = time(step) + dt;
end
% calculate analytical solution below
t=0:0.05:5;
N_analytical=N_nuclei_initial*exp(-t/tau);
% Plot both numerical and analytical solution
plot(time,N_nuclei,'ro',t,N_analytical,'b'); %plots the numerical solution in red and the analytical solution in blue
xlabel('Time (s)')
ylabel('Number of nuclei')
text(2,80,'Time constant = 1s')
text(2,70,'Time step = 0.05s')

運行程序,得到:


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM