1-1 基于matlab的平面曲線曲率的數值計算

1. 工具

   向量函數：設曲線 r ( s ) = ( x ( s ) , y ( s ) ) r(s)=(x(s),y(s)) r(s)=(x(s),y(s))是一條正則曲線，其中 s s s是弧長參數。 r ( s ) r(s) r(s)是以向量形式表示的，是以稱為向量函數。

2. 曲率公式

   • 給定函數為向量函數 r ( s ) = ( x ( s ) , y ( s ) ) r(s)=(x(s),y(s)) r(s)=(x(s),y(s))形式：

     曲率的計算公式為： k = ∣ r ′ ′ ( s ) ∣ (1) k=|r''(s)| \tag {1} k=∣r′′(s)∣(1) 注意 : 當 給定函數 r ( t ) = ( x ( t ) , y ( t ) ) r(t)=(x(t),y(t)) r(t)=(x(t),y(t))參數不是弧長參數時我們需要計算函數的弧長參數，才能進行下一步計算。其中弧長參數計算方法為計算弧微分 s = ∫ 0 t ∣ r ′ ′ ( u ) ∣ d u s=\displaystyle \int^{t}_{0}{|r''(u)|du} s=∫0t​∣r′′(u)∣du.

   • 給定函數為 y = f ( x ) y=f(x) y=f(x)形式：

     曲率 k = ∣ f ′ ′ ( x ) ∣ ( 1 + f ′ 2 ) ( 3 2 ) (2) k=\frac{|f''(x)|}{\sqrt{(1+f'^2)}^{(\frac{3}{2})}}\tag {2} k=(1+f′2)

     ​(23​)∣f′′(x)∣​(2)

3. matlab程式實作

   1.計算圓 r ( t ) = ( 2 c o s t , 2 s i n t ) ， t ∈ ( 0 , 2 π ) r(t)=(2cost,2sint)，t\in (0,2\pi) r(t)=(2cost,2sint)，t∈(0,2π)的曲率

   解: 由公式 ( 1 ) (1) (1)計算曲率時發現參數不是弧長參數，是以我們要計算弧微分 s = ∫ 0 t ∣ r ′ ′ ( u ) ∣ d u s=\displaystyle \int^{t}_{0}{|r''(u)|du} s=∫0t​∣r′′(u)∣du

   于是 s = ∫ 0 t 2 d u = 2 t s=\displaystyle \int^{t}_{0}{2du}=2t s=∫0t​2du=2t

   進而轉換成計算 r ( t ) = ( 2 c o s s 2 , 2 s i n s 2 ) ， s ∈ ( 0 , 4 π ) 的 曲 率 r(t)=(2cos\frac{s}{2},2sin\frac{s}{2})，s\in (0,4\pi)的曲率 r(t)=(2cos2s​,2sin2s​)，s∈(0,4π)的曲率

   MATLAB程式：

clc,clear
h=0.01;                        %步長
s=0:h:4*pi;
x=2*cos(s./2);y=2*sin(s./2);   %定義x，y
r=[x,y];                       %定義向量函數
r1=diff(r)./h;                 %一階導
r2=diff(r1)./h;                %二階導
k=r2;
for i=1:length(r2)
    k(i)=norm(r2(i));
end
           

1. 計算 y = s i n x , x ∈ ( 0 , π ) y=sinx,x\in(0,\pi) y=sinx,x∈(0,π)的曲率

   由公式 ( 2 ) (2) (2) k = ∣ s i n x ∣ ( 1 + c o s x 2 ) ( 3 2 ) k=\frac{|sinx|}{\sqrt{(1+cosx^2)}^{(\frac{3}{2})}} k=(1+cosx2)

   ​(23​)∣sinx∣​

   MATLAB程式：

clc,clear
h=0.01;                        %步長
x=0:h:pi;
y=sin(x);                      %定義x，y
y1=diff(y)./h;                 %一階導
y2=diff(y1)./h;                %二階導
y2(length(y1))=y2(end);        %二階導右端點用向後差商代替，否則會有次元不一緻情況
k=abs(y2)./(sqrt(1+y1.^2).^(3/2));