• Matlab数值计算最简单的一个例子——指数衰减


    放射性衰变是指数衰减的典型例子。另外还有化学反应某反应物的减少,RC电路电流的减小,大气压随海拔高度的减小等。

    指数衰减的方程:

    egin{equation}
    frac{dN(t)}{dt}=-frac{N(t)}{ au}
    label{eq1}
    end{equation}

    其中,(N(t))(t)时刻的物理量(N),对于放射性衰变,(N)就是未衰变的原子核数目。( au)为时间常数。

    方程eqref{eq1}有解析解:

    [N(t)=N(0)exp(-t/ au) ]

    这个解可以通过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}approxfrac{N(t)-N(t-Delta t)}{Delta t} ]

    将以上两式带入方程eqref{eq1},得离散之后的方程:

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

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

    下面写个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')
    

    运行程序,得到:

  • 相关阅读:
    PHP的CURL方法curl_setopt()函数案例介绍(抓取网页,POST数据)
    PHP中CURL方法curl_setopt()函数的参数
    elasticsearch中如何高效的使用filter
    elasticsearch与mongodb分布式集群环境下数据同步
    最完整的Elasticsearch 基础教程
    PHP 解析 ElasticSearch 的 json 方法,有關遍歷所有 json 元素
    分布式搜索引擎Elasticsearch PHP类封装 使用原生api
    [转] Form 表单数据处理 简单教程 formidable 使用心得
    [转] node升级到8.0.0在vscode启动js执行文件报错
    [转] babel的使用
  • 原文地址:https://www.cnblogs.com/joyfulphysics/p/4726503.html
Copyright © 2020-2023  润新知