微分方程

微分方程初步

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

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

Chapter 1: 簡介

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

How to use Euler method to solve the DE一階微分方程，y' = y+hF(x, y)，( Cpp code )， Python code to analyze (py_code)，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 = -2r*y-wwx;

}

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 + hdxdt; // Euler method

ky = y + hdydt;

*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?

(2階微分方程 movie ，Python code )

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

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

dy/dt = -2ry-wwx; // 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，那麼會發生什麼事呢? 讀者是不是開始能夠猜測出振盪變高了，能夠預測出圖形嗎?這樣讀者對微分方程學習有感覺了嗎?

(2階微分方程推導 movie ，Python code )

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] = hdxdt; ry[0] = hdydt; rz[0] = h*dzdt;

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

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

rx[1] = hdxdt; ry[1] = hdydt; rz[1] = h*dzdt;

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

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

rx[2] = hdxdt; ry[2] = hdydt; rz[2] = h*dzdt;

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

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

rx[3] = hdxdt; ry[3] = hdydt; 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)

「彈簧阻尼振盪-微分方程」+ Python code 程式模擬影片(youtube)(pyode)

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

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

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

微分方程

微分方程初步

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

Chapter 1: 簡介

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

How to use Euler method to solve the DE一階微分方程，y' = y+hF(x, y)，( Cpp code )， Python code to analyze (py_code)，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=振盪頻率

那如果當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()

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

dydt = -2r*y-wwx;

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

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

kx = x + hdxdt; // Euler method

ky = y + hdydt;

fclose(pFile); //關檔結束

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

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

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

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)

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

rx[0] = hdxdt; ry[0] = hdydt; rz[0] = h*dzdt;

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

rx[1] = hdxdt; ry[1] = hdydt; rz[1] = h*dzdt;

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

rx[2] = hdxdt; ry[2] = hdydt; rz[2] = h*dzdt;

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

rx[3] = hdxdt; ry[3] = hdydt; 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;