微分方程

 微分方程初步

今天給大家介紹物理學中最常用到的微分方程(Differential Equations)，目的是給有興趣的同學，能迅速掌握物理上微分方程意義，又能避開複雜的數學定義，讓也不是念物理的人知道數學底蘊其實是所有動力學的基礎，認真來說，其實微分方程(DE)真的是一個很有趣的學科，本文編排以微分方程概論為藍圖 中文本 ，《微分方程概論》 林昭仁著。

我的讀書教學筆記，巴斯隨筆， 20200410

Chapter 1: 簡介  

Chapter 2: 一階常微分方程 (ref)

<< Cpp code >>How to use Euler method to solve the DE一階微分方程，y' = y+hF(x, y) ，AI LSTM  jupyter code (lstm_code)

// Eulers Method to solve a differential equation

#include<iostream>  

#include<iomanip>

#include<cmath>       // math 

#define hmt 1000      // run steps

using namespace std;

void spring( double x, double y, double *dxdt,double *dydt)  // spring DEs

{

    const double r=0.0, w=2.0;

   *dxdt = y;

   *dydt = -2*r*y-w*w*x;

}


void Euler(double h, void(*f)( double x, double y, double *dxdt,double *dydt),double *x,double *y)  // 函式指標 

{

   double kx, ky;

   double dxdt, dydt;

   f(*x,*y,&dxdt,&dydt);

   kx = *x + h*dxdt;              // Euler method

   ky = *y + h*dydt;

   *x = kx;

   *y = ky;

}


void main()   // 主程式

{

   double h = 0.01;  // h間隔，h 太大，誤差變大

   FILE *pFile;          // 開檔指標

   pFile = fopen("r_0_w_2.txt","w");  // 指定存檔檔名

   int i;

   double x = 1, y = 0;

   for(i=0; i < hmt; i++)

   {

      fprintf(pFile, "%10d %12.6f %12.6f\n",i,x,y);  // 寫檔

      Euler(h, spring, &x, &y);     // 呼叫Euler method

   }

   fclose(pFile);  //關檔  結束

}

 Question 1. 彈簧拉動物體運動的問題，F=ma，若恢復力為彈簧彈力，F=-kx，k>0，再考慮阻尼力，f=-bv，b為參數，v為速度，如何列出ODE?

如果用Euler method來解微分方程ODE，運動微分程式如上面所示，利用變數變換，二階微分程式拆成兩個一階微分程式組，再用Euler method數值解此兩個一階微分程式組，分析如下: 

 dx/dt = y;                           // x 是位移，y 是速度

 dy/dt = -2ry-w*w*x;  // b= 2r=阻尼，w=振盪頻率

如上圖，是利用此程式產生數據，excel畫出，r=0視為無阻尼，視為標準簡諧振盪，當 h= 0.0001 間隔比較小，y 對 x 為正圓振動相圖，x, y 相位相差90度，但當間隔 h= 0.01， 會發生什麼事?這裡的x是位移，y代表的是速度項，二階微分方程，改成兩個一階微分方程組來解，這是個小技巧，如程式所示。

間隔 h = 0.01，r=0視為無阻尼，應該為標準簡諧振盪，但是y 對 x 為正圓振動相圖無法密合，這是使用Euler method 常見的缺點之一，改進之道為使用runge Kutta 4th 方法。或是h縮小，那如果當b=2r=2， 振盪頻率 w= 4，那麼會發生什麼事呢? 建議讀者，設定各種參數(r, w)，產生各種數值畫圖，這樣對微分方程學習與掌握會更有感覺。

那如果當b=2r=2， 振盪頻率 w= 10，那麼會發生什麼事呢? 讀者是不是開始能夠猜測出振盪變高了，能夠預測出圖形嗎?這樣讀者對微分方程學習有感覺了嗎?

Lorenz code by using Runge Kutta 4th (Lorenz.cpp  written by PHLee)

#include<stdio.h>

#define hmt 2000

void lorenz( double t, double x, double y, double z, double *dxdt, double *dydt

 , double *dzdt)

{

      const double sigma=10.0, b=8.0/3.0, r=28.0;

     *dxdt = -sigma * x + sigma * y;

    *dydt = -x * z +r * x - y;

    *dzdt = x * y - b * z;

}


void r_k(double h, double t, void(*f)( double t, double x, double y, double z, double *dxdt, double *dydt  , double *dzdt)

 ,double *x, double *y, double *z)

 {

    double kx, ky, kz;

    double dxdt, dydt, dzdt;

    double rx[4], ry[4], rz[4];

    f(t, *x, *y, *z, &dxdt, &dydt, &dzdt);

    rx[0] = h*dxdt; ry[0] = h*dydt; rz[0] = h*dzdt;

    kx = *x+0.5*rx[0]; ky = *y+0.5*ry[0]; kz = *z+0.5*rz[0];


    f(t+0.5*h, kx, ky, kz, &dxdt, &dydt, &dzdt);

    rx[1] = h*dxdt; ry[1] = h*dydt; rz[1] = h*dzdt;

    kx = *x+0.5*rx[1]; ky = *y+0.5*ry[1]; kz = *z+0.5*rz[1];


    f(t+0.5*h, kx, ky, kz, &dxdt, &dydt, &dzdt);

    rx[2] = h*dxdt; ry[2] = h*dydt; rz[2] = h*dzdt;

    kx = *x+0.5*rx[2]; ky = *y+0.5*ry[2]; kz = *z+0.5*rz[2];


    f(t+h, kx, ky, kz, &dxdt, &dydt, &dzdt);

    rx[3] = h*dxdt; ry[3] = h*dydt; rz[3] = h*dzdt;

    *x = *x + (rx[0] + rx[1]*2 + rx[2]*2 + rx[3])/6.0;

    *y = *y + (ry[0] + ry[1]*2 + ry[2]*2 + ry[3])/6.0;

    *z = *z + (rz[0] + rz[1]*2 + rz[2]*2 + rz[3])/6.0;

}


 int main()

 {

    double h=0.01;

    FILE *pFile;

    pFile = fopen ("h_2000_s_10_r_28_b_8_3/s_10_r_28_b_8_3.txt", "w") ;

    int n, t=1, i;

    double vx = 10.0, vy = 20.0, vz = 30.0;

    for(i=0; i<hmt; i++)

    {

       fprintf(pFile,"%f,%f,%f\n",vx, vy, vz);

       r_k(h, h*t, lorenz, &vx, &vy, &vz);

   }

  fclose(pFile);

  return 0;

  }




微積分入門(李柏翰老師教學影片) 

微分方程教學影片(youtube)

20210713 物理培訓 (ppt1)(ppt2)(微分方程PPT教學影片)

20210806  中原大學資工系 黃琮暐教授

主題： 量子電腦的優勢 (movie)