Embedded Files
微分方程初步
微分方程初步
今天給大家介紹物理學中最常用到的微分方程(Differential Equations)，目的是給有興趣的同學，能迅速掌握物理上微分方程意義，又能避開複雜的數學定義，讓也不是念物理的人知道數學底蘊其實是所有動力學的基礎，認真來說，其實微分方程(DE)真的是一個很有趣的學科，本文編排以微分方程概論為藍圖 中文本 ，《微分方程概論》 林昭仁著。
今天給大家介紹物理學中最常用到的微分方程(Differential Equations)，目的是給有興趣的同學，能迅速掌握物理上微分方程意義，又能避開複雜的數學定義，讓也不是念物理的人知道數學底蘊其實是所有動力學的基礎，認真來說，其實微分方程(DE)真的是一個很有趣的學科，本文編排以微分方程概論為藍圖 中文本 ，《微分方程概論》 林昭仁著。
我的讀書教學筆記，巴斯隨筆， 20200410
我的讀書教學筆記，巴斯隨筆， 20200410
Chapter 1: 簡介
Chapter 1: 簡介
How to use Euler method to solve the DE一階微分方程，y' = y+hF(x, y)，( Cpp code )， Python code to analyze (py_code)(py_code_v2)，AI LSTM jupyter code (lstm_code)
How to use Euler method to solve the DE一階微分方程，y' = y+hF(x, y)，( Cpp code )， Python code to analyze (py_code)(py_code_v2)，AI LSTM jupyter code (lstm_code)
// Eulers Method to solve a differential equation
// Eulers Method to solve a differential equation
#include<iostream>
#include<iostream>
#include<iomanip>
#include<iomanip>
#include<cmath> // math
#include<cmath> // math
#define hmt 1000 // run steps
#define hmt 1000 // run steps
using namespace std;
using namespace std;
void spring( double x, double y, double *dxdt,double *dydt) // spring DEs
void spring( double x, double y, double *dxdt,double *dydt) // spring DEs
{
{
const double r=0.0, w=2.0;
const double r=0.0, w=2.0;
*dxdt = y;
*dxdt = y;
*dydt = -2*r*y-w*w*x;
*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) // 函式指標
void Euler(double h, void(*f)( double x, double y, double *dxdt,double *dydt),double *x,double *y) // 函式指標
{
{
double kx, ky;
double kx, ky;
double dxdt, dydt;
double dxdt, dydt;
f(*x,*y,&dxdt,&dydt);
f(*x,*y,&dxdt,&dydt);
kx = *x + h*dxdt; // Euler method
kx = *x + h*dxdt; // Euler method
ky = *y + h*dydt;
ky = *y + h*dydt;
*x = kx;
*x = kx;
*y = ky;
*y = ky;
}
}
void main() // 主程式
void main() // 主程式
{
{
double h = 0.01; // h間隔，h 太大，誤差變大
double h = 0.01; // h間隔，h 太大，誤差變大
FILE *pFile; // 開檔指標
FILE *pFile; // 開檔指標
pFile = fopen("r_0_w_2.txt","w"); // 指定存檔檔名
pFile = fopen("r_0_w_2.txt","w"); // 指定存檔檔名
int i;
int i;
double x = 1, y = 0;
double x = 1, y = 0;
for(i=0; i < hmt; i++)
for(i=0; i < hmt; i++)
{
{
fprintf(pFile, "%10d %12.6f %12.6f\n",i,x,y); // 寫檔
fprintf(pFile, "%10d %12.6f %12.6f\n",i,x,y); // 寫檔
Euler(h, spring, &x, &y); // 呼叫Euler method
Euler(h, spring, &x, &y); // 呼叫Euler method
}
}
fclose(pFile); //關檔 結束
fclose(pFile); //關檔 結束
}
}
Question 1. 彈簧拉動物體運動的問題，F=ma，若恢復力為彈簧彈力，F=-kx，k>0，再考慮阻尼力，f=-bv，b為參數，v為速度，如何列出ODE?
Question 1. 彈簧拉動物體運動的問題，F=ma，若恢復力為彈簧彈力，F=-kx，k>0，再考慮阻尼力，f=-bv，b為參數，v為速度，如何列出ODE?
(2階微分方程 movie ，Python code )
如果用Euler method來解微分方程ODE，運動微分程式如上面所示，利用變數變換，二階微分程式拆成兩個一階微分程式組，再用Euler method數值解此兩個一階微分程式組，分析如下:
如果用Euler method來解微分方程ODE，運動微分程式如上面所示，利用變數變換，二階微分程式拆成兩個一階微分程式組，再用Euler method數值解此兩個一階微分程式組，分析如下:
dx/dt = y; // x 是位移，y 是速度
dx/dt = y; // x 是位移，y 是速度
dy/dt = -2ry-w*w*x; // b= 2r=阻尼，w=振盪頻率
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代表的是速度項，二階微分方程，改成兩個一階微分方程組來解，這是個小技巧，如程式所示。
如上圖，是利用此程式產生數據，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)，產生各種數值畫圖，這樣對微分方程學習與掌握會更有感覺。
間隔 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，那麼會發生什麼事呢? 讀者是不是開始能夠猜測出振盪變高了，能夠預測出圖形嗎?這樣讀者對微分方程學習有感覺了嗎?
那如果當b=2r=2， 振盪頻率 w= 10，那麼會發生什麼事呢? 讀者是不是開始能夠猜測出振盪變高了，能夠預測出圖形嗎?這樣讀者對微分方程學習有感覺了嗎?
(2階微分方程推導 movie ，Python code )
#include<stdio.h>
#include<stdio.h>
#define hmt 2000
#define hmt 2000
void lorenz( double t, double x, double y, double z, double *dxdt, double *dydt
void lorenz( double t, double x, double y, double z, double *dxdt, double *dydt
, double *dzdt)
, double *dzdt)
{
{
const double sigma=10.0, b=8.0/3.0, r=28.0;
const double sigma=10.0, b=8.0/3.0, r=28.0;
*dxdt = -sigma * x + sigma * y;
*dxdt = -sigma * x + sigma * y;
*dydt = -x * z +r * x - y;
*dydt = -x * z +r * x - y;
*dzdt = x * y - b * z;
*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)
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 *x, double *y, double *z)
{
{
double kx, ky, kz;
double kx, ky, kz;
double dxdt, dydt, dzdt;
double dxdt, dydt, dzdt;
double rx[4], ry[4], rz[4];
double rx[4], ry[4], rz[4];
f(t, *x, *y, *z, &dxdt, &dydt, &dzdt);
f(t, *x, *y, *z, &dxdt, &dydt, &dzdt);
rx[0] = h*dxdt; ry[0] = h*dydt; rz[0] = h*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];
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);
f(t+0.5*h, kx, ky, kz, &dxdt, &dydt, &dzdt);
rx[1] = h*dxdt; ry[1] = h*dydt; rz[1] = h*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];
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);
f(t+0.5*h, kx, ky, kz, &dxdt, &dydt, &dzdt);
rx[2] = h*dxdt; ry[2] = h*dydt; rz[2] = h*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];
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);
f(t+h, kx, ky, kz, &dxdt, &dydt, &dzdt);
rx[3] = h*dxdt; ry[3] = h*dydt; rz[3] = h*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;
*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;
*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;
*z = *z + (rz[0] + rz[1]*2 + rz[2]*2 + rz[3])/6.0;
}
}
int main()
int main()
{
{
double h=0.01;
double h=0.01;
FILE *pFile;
FILE *pFile;
pFile = fopen ("h_2000_s_10_r_28_b_8_3/s_10_r_28_b_8_3.txt", "w") ;
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;
int n, t=1, i;
double vx = 10.0, vy = 20.0, vz = 30.0;
double vx = 10.0, vy = 20.0, vz = 30.0;
for(i=0; i<hmt; i++)
for(i=0; i<hmt; i++)
{
{
fprintf(pFile,"%f,%f,%f\n",vx, vy, vz);
fprintf(pFile,"%f,%f,%f\n",vx, vy, vz);
r_k(h, h*t, lorenz, &vx, &vy, &vz);
r_k(h, h*t, lorenz, &vx, &vy, &vz);
}
}
fclose(pFile);
fclose(pFile);
return 0;
return 0;
}
}
Page updated
Google Sites
Report abuse