MATLABのdiffコマンドでシンボリック微分を計算

MATLABには、シンボリック微分を計算するdiffコマンドがあります。最も簡単な形式では、微分したい関数をdiffコマンドの引数として渡します。

例えば、関数f(t) = 3t² + 2t^-2の導関数を計算してみましょう。

例

スクリプトファイルを作成し、以下のコードを入力します。

 syms t f = 3*t^2 + 2*t^(-2); diff(f) 

上記のコードをコンパイルして実行すると、以下の結果が得られます。

 ans = 6*t - 4/t^3 

以下は上記計算のOctave版です。

 pkg load symbolic symbols t = sym("t"); f = 3*t^2 + 2*t^(-2); differentiate(f,t) 

Octaveはコードを実行し、以下の結果を返します。

 ans = -(4.0)*t^(-3.0)+(6.0)*t 

微分の基本ルールの検証

関数の微分に関するさまざまな方程式やルールを簡単に述べ、それらを検証します。この目的で、1階導関数をf'(x)、2階導関数をf"(x)で表します。

以下は微分のルールです。

ルール1

任意の関数f、gと任意の実数a、bに対して、関数

h(x) = af(x) + bg(x)のxに関する導関数は

h'(x) = af'(x) + bg'(x)

ルール2

和と差のルールでは、fとgが2つの関数、f'とg'がそれぞれの導関数のとき、

(f + g)' = f' + g'

(f - g)' = f' - g'

ルール3

積のルールでは、fとgが2つの関数、f'とg'がそれぞれの導関数のとき、

(f·g)' = f'·g + g'·f

ルール4

商のルールでは、fとgが2つの関数、f'とg'がそれぞれの導関数のとき、

(f/g)' = (f'·g - g'·f)/g²

ルール5

多項式または基本的な冪乗ルールでは、y = f(x) = xⁿのとき、f' = n·x^(n-1)

このルールの直接的な帰結として、任意の定数の導関数はゼロです。つまり、y = k（任意の定数）のとき、

f' = 0

ルール6

連鎖ルールでは、合成関数h(x) = f(g(x))のxに関する導関数は、

h'(x) = f'(g(x))·g'(x)

例

スクリプトファイルを作成し、以下のコードを入力します。

 syms x syms t f = (x + 2)*(x^2 + 3) der1 = diff(f) f = (t^2 + 3)*(sqrt(t) + t^3) der2 = diff(f) f = (x^2 - 2*x + 1)*(3*x^3 - 5*x^2 + 2) der3 = diff(f) f = (2*x^2 + 3*x)/(x^3 + 1) der4 = diff(f) f = (x^2 + 1)^17 der5 = diff(f) f = (t^3 + 3* t^2 + 5*t -9)^(-6) der6 = diff(f) 

ファイルを実行すると、MATLABは以下の結果を表示します。

 f = (x^2 + 3)*(x + 2) der1 = 2*x*(x + 2) + x^2 + 3 f = (t^(1/2) + t^3)*(t^2 + 3) der2 = (t^2 + 3)*(3*t^2 + 1/(2*t^(1/2))) + 2*t*(t^(1/2) + t^3) f = (x^2 - 2*x + 1)*(3*x^3 - 5*x^2 + 2) der3 = (2*x - 2)*(3*x^3 - 5*x^2 + 2) - (- 9*x^2 + 10*x)*(x^2 - 2*x + 1) f = (2*x^2 + 3*x)/(x^3 + 1) der4 = (4*x + 3)/(x^3 + 1) - (3*x^2*(2*x^2 + 3*x))/(x^3 + 1)^2 f = (x^2 + 1)^17 der5 = 34*x*(x^2 + 1)^16 f = 1/(t^3 + 3*t^2 + 5*t - 9)^6 der6 = -(6*(3*t^2 + 6*t + 5))/(t^3 + 3*t^2 + 5*t - 9)^7 

以下は上記計算のOctave版です。

 pkg load symbolic symbols x = sym("x"); t = sym("t"); f = (x + 2)*(x^2 + 3) der1 = differentiate(f,x) f = (t^2 + 3)*(t^(1/2) + t^3) der2 = differentiate(f,t) f = (x^2 - 2*x + 1)*(3*x^3 - 5*x^2 + 2) der3 = differentiate(f,x) f = (2*x^2 + 3*x)/(x^3 + 1) der4 = differentiate(f,x) f = (x^2 + 1)^17 der5 = differentiate(f,x) f = (t^3 + 3* t^2 + 5*t -9)^(-6) der6 = differentiate(f,t) 

Octaveはコードを実行し、以下の結果を返します。

 f = (2.0+x)*(3.0+x^(2.0)) der1 = 3.0+x^(2.0)+(2.0)*(2.0+x)*x f = (t^(3.0)+sqrt(t))*(3.0+t^(2.0)) der2 = (2.0)*(t^(3.0)+sqrt(t))*t+((3.0)*t^(2.0)+(0.5)*t^(-0.5))*(3.0+t^(2.0)) f = (1.0+x^(2.0)-(2.0)*x)*(2.0-(5.0)*x^(2.0)+(3.0)*x^(3.0)) der3 = (-2.0+(2.0)*x)*(2.0-(5.0)*x^(2.0)+(3.0)*x^(3.0))+((9.0)*x^(2.0)-(10.0)*x)*(1.0+x^(2.0)-(2.0)*x) f = (1.0+x^(3.0))^(-1)*((2.0)*x^(2.0)+(3.0)*x) der4 = (1.0+x^(3.0))^(-1)*(3.0+(4.0)*x)-(3.0)*(1.0+x^(3.0))^(-2)*x^(2.0)*((2.0)*x^(2.0)+(3.0)*x) f = (1.0+x^(2.0))^(17.0) der5 = (34.0)*(1.0+x^(2.0))^(16.0)*x f = (-9.0+(3.0)*t^(2.0)+t^(3.0)+(5.0)*t)^(-6.0) der6 = -(6.0)*(-9.0+(3.0)*t^(2.0)+t^(3.0)+(5.0)*t)^(-7.0)*(5.0+(3.0)*t^(2.0)+(6.0)*t) 

指数関数、対数関数、三角関数の導関数

以下の表は、よく使用される指数関数、対数関数、三角関数の導関数を示しています。

関数	導関数
ca·x	ca·x·ln c·a（lnは自然対数）
ex	ex
ln x	1/x
logcx	1/(x·ln c)
xx	xx·(1 + ln x)
sin(x)	cos(x)
cos(x)	-sin(x)
tan(x)	sec2(x)、または1/cos2(x)、または1 + tan2(x)
cot(x)	-csc2(x)、または-1/sin2(x)、または-(1 + cot2(x))
sec(x)	sec(x)·tan(x)
csc(x)	-csc(x)·cot(x)

例

スクリプトファイルを作成し、以下のコードを入力します。

 syms x y = exp(x) diff(y) y = x^9 diff(y) y = sin(x) diff(y) y = tan(x) diff(y) y = cos(x) diff(y) y = log(x) diff(y) y = log10(x) diff(y) y = sin(x)^2 diff(y) y = cos(3*x^2 + 2*x + 1) diff(y) y = exp(x)/sin(x) diff(y) 

ファイルを実行すると、MATLABは以下の結果を表示します。

 y = exp(x) ans = exp(x) y = x^9 ans = 9*x^8 y = sin(x) ans = cos(x) y = tan(x) ans = tan(x)^2 + 1 y = cos(x) ans = -sin(x) y = log(x) ans = 1/x y = log(x)/log(10) ans = 1/(x*log(10)) y = sin(x)^2 ans = 2*cos(x)*sin(x) y = cos(3*x^2 + 2*x + 1) ans = -sin(3*x^2 + 2*x + 1)*(6*x + 2) y = exp(x)/sin(x) ans = exp(x)/sin(x) - (exp(x)*cos(x))/sin(x)^2 

以下は上記計算のOctave版です。

 pkg load symbolic symbols x = sym("x"); y = Exp(x) differentiate(y,x) y = x^9 differentiate(y,x) y = Sin(x) differentiate(y,x) y = Tan(x) differentiate(y,x) y = Cos(x) differentiate(y,x) y = Log(x) differentiate(y,x) % symbolicパッケージはこれ未対応 %y = Log10(x) %differentiate(y,x) y = Sin(x)^2 differentiate(y,x) y = Cos(3*x^2 + 2*x + 1) differentiate(y,x) y = Exp(x)/Sin(x) differentiate(y,x) 

Octaveはコードを実行し、以下の結果を返します。

 y = exp(x) ans = exp(x) y = x^(9.0) ans = (9.0)*x^(8.0) y = sin(x) ans = cos(x) y = tan(x) ans = 1+tan(x)^2 y = cos(x) ans = -sin(x) y = log(x) ans = x^(-1) y = sin(x)^(2.0) ans = (2.0)*sin(x)*cos(x) y = cos(1.0+(2.0)*x+(3.0)*x^(2.0)) ans = -(2.0+(6.0)*x)*sin(1.0+(2.0)*x+(3.0)*x^(2.0)) y = sin(x)^(-1)*exp(x) ans = sin(x)^(-1)*exp(x)-sin(x)^(-2)*cos(x)*exp(x) 

高階導関数の計算

関数fの高階導関数を計算するには、構文diff(f,n)を使用します。

関数y = f(x) = x·e^-3xの2階導関数を計算してみましょう。

 f = x*exp(-3*x); diff(f, 2) 

MATLABはコードを実行し、以下の結果を返します。

 ans = 9*x*exp(-3*x) - 6*exp(-3*x) 

以下は上記計算のOctave版です。

 pkg load symbolic symbols x = sym("x"); f = x*Exp(-3*x); differentiate(f, x, 2) 

Octaveはコードを実行し、以下の結果を返します。

 ans = (9.0)*exp(-(3.0)*x)*x-(6.0)*exp(-(3.0)*x) 

例

この例では、ある問題を解きます。関数y = f(x) = 3 sin(x) + 7 cos(5x)が与えられています。方程式f" + f = -5cos(2x)が成り立つかどうかを調べます。

スクリプトファイルを作成し、以下のコードを入力します。

 syms x y = 3*sin(x)+7*cos(5*x); % 関数を定義 lhs = diff(y,2)+y; % 方程式の左辺を評価 rhs = -5*cos(2*x); % 方程式の右辺 if(isequal(lhs,rhs)) disp('Yes, the equation holds true'); else disp('No, the equation does not hold true'); end disp('Value of LHS is: '), disp(lhs); 

ファイルを実行すると、以下の結果が表示されます。

 No, the equation does not hold true Value of LHS is: -168*cos(5*x) 

以下は上記計算のOctave版です。

 pkg load symbolic symbols x = sym("x"); y = 3*Sin(x)+7*Cos(5*x); % 関数を定義 lhs = differentiate(y, x, 2) + y; % 方程式の左辺を評価 rhs = -5*Cos(2*x); % 方程式の右辺 if(lhs == rhs) disp('Yes, the equation holds true'); else disp('No, the equation does not hold true'); end disp('Value of LHS is: '), disp(lhs); 

Octaveはコードを実行し、以下の結果を返します。

 No, the equation does not hold true Value of LHS is: -(168.0)*cos((5.0)*x) 

曲線の極大値と極小値の発見

グラフの局所的な極大値と極小値を探す場合、基本的に関数のグラフ上で特定の局所、またはシンボリック変数の特定の値の範囲での最高点や最低点を探しています。

関数y = f(x)において、グラフの傾きがゼロになる点を定常点と呼びます。言い換えると、定常点はf'(x) = 0となる点です。

関数の定常点を見つけるには、微分して導関数をゼロに等しく設定し、方程式を解きます。

例

関数f(x) = 2x³ + 3x² − 12x + 17の定常点を求めます。

以下の手順を実行します。

まず関数を入力し、そのグラフをプロットします。

 syms x y = 2*x^3 + 3*x^2 - 12*x + 17; % 関数を定義 ezplot(y) 

MATLABはコードを実行し、以下のプロットを返します。

上記例のOctave版コードは以下の通りです。

 pkg load symbolic symbols x = sym('x'); y = inline("2*x^3 + 3*x^2 - 12*x + 17"); ezplot(y) print -deps graph.eps 

グラフ上で局所的な極大値と極小値を見つけることが目的なので、グラフ上の区間[-2, 2]で局所的な極大値と極小値を求めます。

 syms x y = 2*x^3 + 3*x^2 - 12*x + 17; % 関数を定義 ezplot(y, [-2, 2]) 

MATLABはコードを実行し、以下のプロットを返します。

上記例のOctave版コードは以下の通りです。

 pkg load symbolic symbols x = sym('x'); y = inline("2*x^3 + 3*x^2 - 12*x + 17"); ezplot(y, [-2, 2]) print -deps graph.eps 

次に導関数を計算します。

 g = diff(y) 

MATLABはコードを実行し、以下の結果を返します。

 g = 6*x^2 + 6*x - 12 

上記計算のOctave版は以下の通りです。

 pkg load symbolic symbols x = sym("x"); y = 2*x^3 + 3*x^2 - 12*x + 17; g = differentiate(y,x) 

Octaveはコードを実行し、以下の結果を返します。

 g = -12.0+(6.0)*x+(6.0)*x^(2.0) 

導関数関数gを解いて、ゼロになる値を求めます。

 s = solve(g) 

MATLABはコードを実行し、以下の結果を返します。

 s = 1 -2 

上記計算のOctave版は以下の通りです。

 pkg load symbolic symbols x = sym("x"); y = 2*x^3 + 3*x^2 - 12*x + 17; g = differentiate(y,x) roots([6, 6, -12]) 

Octaveはコードを実行し、以下の結果を返します。

 g = -12.0+(6.0)*x^(2.0)+(6.0)*x ans = -2 1 

これはプロットと一致します。では、臨界点x = 1、-2で関数fを評価します。 シンボリック関数に値を代入するには、subsコマンドを使用します。

 subs(y, 1), subs(y, -2) 

MATLABはコードを実行し、以下の結果を返します。

 ans = 10 ans = 37 

上記計算のOctave版は以下の通りです。

 pkg load symbolic symbols x = sym("x"); y = 2*x^3 + 3*x^2 - 12*x + 17; g = differentiate(y,x) roots([6, 6, -12]) subs(y, x, 1), subs(y, x, -2) 

 ans = 10.0 ans = 37.0-4.6734207789940138748E-18*I 

したがって、区間[-2,2]における関数f(x) = 2x³ + 3x² − 12x + 17の最小値と最大値は、それぞれ10と37です。

微分方程式の解法

MATLABには、微分方程式をシンボリックに解くdsolveコマンドがあります。

単一の方程式の解を見つけるためのdsolveコマンドの最も基本的な形式は

 dsolve('eqn') 

ここで、eqnは方程式を入力するためのテキスト文字列です。

MATLABがC1、C2などとラベル付けする任意定数のセットを含むシンボリック解を返します。

問題の初期条件と境界条件を、方程式に続くカンマ区切りのリストとして指定することもできます。

 dsolve('eqn','cond1', 'cond2',…) 

dsolveコマンドを使用する目的で、導関数はDで示します。例えば、f'(t) = -2*f + cos(t)のような方程式は以下のように入力します。

'Df = -2*f + cos(t)'

高階導関数は、Dに続く導関数の次数で示します。

例えば、f"(x) + 2f'(x) = 5sin3xの方程式は以下のように入力します。

'D2y + 2Dy = 5*sin(3*x)'

1階微分方程式の簡単な例：y' = 5yを取り上げます。

 s = dsolve('Dy = 5*y') 

MATLABはコードを実行し、以下の結果を返します。

 s = C2*exp(5*t) 

2階微分方程式の別の例：y" - y = 0、y(0) = -1、y'(0) = 2。

 dsolve('D2y - y = 0','y(0) = -1','Dy(0) = 2') 

MATLABはコードを実行し、以下の結果を返します。

 ans = exp(t)/2 - (3*exp(-t))/2