ここに書いたコードはR2021aで実行しています.
はじめに
COVID-19の影響により, 学校でのテスト形態も様変わりしてきました.
今までは「計算ミスをしても途中式を見てくれたからある程度点が取れていた」という人も, 「テストがオンラインで答えのみを入力する方式になってしまったから, 勉強した内容をしっかり理解していても, 点が取れない」という経験をよく耳にします (耳にするだけでなく, ごく最近, 自分もやらかしました..)
この記事は, そんな「理解度90%以上なのに, いい点数が取れない」という方に届けたく, 書き始めました.
私たちはコンピュータではないので, 間違えるのは当たり前です!!
もし, 手計算に自信がない場合は以下の方法で自分の答えを確認してみてはいかがでしょうか!!
シリーズ他記事
1. 線形代数編
MATLABで数式を扱う
今回は, 微分積分の数式をコンピュータ内で扱うために, MATLAB Symbolic math Toolbox を使います.
以下の式を例に, Symbolic Math Toolboxの簡単な使い方について書いていこうと思います.
問題
$x=2$での$y$の値を求めなさい!!
$$ y = 3x^2+2x+1 $$
数式を作る
上の例を作るには, $x$というSymbolic変数を作る必要がある. そのためには, symsを使い, 以下のように変数を宣言する.
syms x
そして, $y$に上の式を代入します.
y = 3*x^2 + 2*x + 1
MATLABではコードの後ろに;をつけなければ,
y =
3*x^2 + 2*x + 1
という風に, その変数の内容を表示してくれます. これで, MATLAB内に数式を作れました.
値の代入
次に, 答えを求めます. 答えの求め方は, $x=2$を代入すれば終わりです. MATLABで数式に値を代入する方法は, 関数subsを用いて,
答え = subs(数式,代入したいSymbolic変数,代入したい値);
とすると, 計算できます. よって,
answer = subs(y,x,2);
answer =
17
となり, 問題が解けました.
(ついでに) 式の可視化
数式を可視化できると, 答えがあっているかだいたいわかってくることはよくあります. 加えて初見の関数の形もだいたい想像できるようになるとか聞いたこともあります. 今回の場合は,
fplot(y)
grid on
だけで,
という感じで, 簡単に作れちゃいます!!
微分積分学
今回は, 理系大学生が必ず通らなければならない科目の1つであろう微分積分学を題材に, 例題を解きながらSymbolic Math Toolboxを使っていきたいと思います!!
極限
微分を語ろうとするうえで避けては通れないのが, 「極限」です.
例題1
問題
以下を求めよ.
$$\lim_{x \rightarrow \infty} x^{\frac{1}{x}}$$
解法
この問題をMATLABで解く場合には,
- 数式$x^{\frac{1}{x}}$ を作成
- 関数limitで極限を求める
答え = limit(数式,操作する変数,求める極限)
という手順になります.
syms x % Symbolic変数xを定義
limex = x^(1/x) % 数式を決定
limit(limex,x,Inf) % 極限を求める
%Visualize (可視化)
figure
fplot(limex)
grid on
xlim([0.1 300]);
limex =
x^(1/x)
ans =
1
例題2
問題
以下を求めよ.
$$\lim_{x \rightarrow 0} \frac{\log (1+x)-\sin x}{x^2}$$
解法
この問題も上の例題1と全く同じ方法で計算できます.
syms x % Symbolic変数xを定義
limex = (log(1+x)-sin(x))/x^2 % 数式を決定
limit(limex,x,0) % 極限を求める
%Visualize
figure
fplot(limex)
grid on
limex =
(log(x + 1) - sin(x))/x^2
ans =
-1/2
微分
さあ, 本題の1つの「微分」です.
１変数関数の微分
問題1
$$f(x) = x^22^x$$
の導関数$\frac{df}{dx}$及び2階の導関数$\frac{d^2f}{dx^2}$を求めよ.
解法1
MATLABでの微分はdiff関数によって,
答え = diff(数式, 微分する変数, 微分の階数)
と求めることができます. これを使うと, 今回の問題は, 以下のように解けます.
syms x
diffex = x^2*2^x
diffex_1 = diff(diffex,x)
diffex_2 = diff(diffex,x,2)
% Visualize
figure
fplot(diffex)
hold on
fplot(diffex_1)
fplot(diffex_2)
hold off
grid on
xlim([-1 1]);
legend('f(x)','df/dx','df^2/dx^2','Location','northwest')
saveas(gcf,'diff1.png');
diffex =
2^x*x^2
diffex_1 =
2*2^x*x + 2^x*x^2*log(2)
diffex_2 =
2*2^x + 4*2^x*x*log(2) + 2^x*x^2*log(2)^2
偏微分
問題2
$$f(x,y) = \frac{x^4-3x^2y^2}{2x^2+y^2}$$
の偏導関数$\frac{\partial f}{\partial x}$と$\frac{\partial f}{\partial y}$を求めよ
解法2
偏微分でも解き方は全く変わりません. また, 2変数を含む数式を可視化するときは, fsurf関数を用いるとよくわかる.
syms x y % Symbolic変数を定義
pdiffex = (x^4-3*x^2*y^2)/(2*x^2+y^2) % 数式を決定
pdiffex_x = diff(pdiffex,x,1) % xによる偏微分を計算
pdiffex_y = diff(pdiffex,y,1) % yによる偏微分を計算
% Visualize
figure
subplot(3,1,1)
fsurf(pdiffex)
grid on
title('f')
subplot(3,1,2)
fsurf(pdiffex_x)
grid on
title('\partial f/\partial x')
subplot(3,1,3)
fsurf(pdiffex_y)
grid on
title('\partial f/\partial y')
pdiffex =
-(- x^4 + 3*x^2*y^2)/(2*x^2 + y^2)
pdiffex_x =
(4*x*(- x^4 + 3*x^2*y^2))/(2*x^2 + y^2)^2 - (- 4*x^3 + 6*x*y^2)/(2*x^2 + y^2)
pdiffex_y =
(2*y*(- x^4 + 3*x^2*y^2))/(2*x^2 + y^2)^2 - (6*x^2*y)/(2*x^2 + y^2)
テーラー展開
微分を使う定理の1つに「テーラー展開」があります. あらゆる近似や数値解析で基礎となる大事な定理です.
問題3
$$f1(x) = \sin(x)$$
$$f2(x) = \cos(x)$$
のテーラー展開を求めよ
解法3
MATLABではテーラー展開をtaylor関数で計算することができる.
syms x % Symbolic変数を定義
f1 = sin(x) % 数式1を決定
f2 = cos(x) % 数式2を決定
T1 = taylor(f1,x,'Order',7) % テーラー展開を7次まで求解
T2 = taylor(f2,x,'Order',9) % テーラー展開を9次まで求解
% Visualize
figure
subplot(2,1,1)
fplot(f1)
hold on
fplot(T1)
hold off
grid on
title('sin(x)')
legend('Original','Taylor');
subplot(2,1,2)
fplot(f2)
hold on
fplot(T2)
hold off
grid on
title('cos(x)')
legend('Original','Taylor');
saveas(gcf,'diff3.png');
f1 =
sin(x)
f2 =
cos(x)
T1 =
x^5/120 - x^3/6 + x
T2 =
x^8/40320 - x^6/720 + x^4/24 - x^2/2 + 1
ちなみに, こんなに簡単にテーラー展開を求められると, 以下のようなこともできちゃいます.
syms x
f = sin(x); % 元の関数
l = 15; % ビデオの長さ
T = taylor(f,x,'Order',1);
v = VideoWriter('taylor.mp4','MPEG-4');
v.FrameRate = 2;
v.Quality = 90;
open(v);
figure
fplot(f);
hold on
pl = fplot(T);
hold off
grid on
ylim([-5 5]);
legend('Original','Taylor');
for i=2:l
T = taylor(f,x,'Order',i);
hold on
pl.Function = T;
str = sprintf('Taylor sin (order=%d)',i);
title(str);
frame = getframe(gcf);
writeVideo(v,frame);
end
close(v);
テーラー展開は次数を増やすほど, より広い範囲を予測してくれます!!
— ヒビヤ@MATLAB Ambassador (@hibs_MATLAB_Amb) June 24, 2021
これが, 常微分方程式数値解法の基礎であり, 有名なルンゲクッタ法もここから生まれたと考えると, 偉大な定理です.
(解くのめちゃくちゃ面倒で, テストには絶対に出てほしくなかったですけど..)#MATLAB #uec #電通大 #数値解析 pic.twitter.com/bq1QGSBxT4
積分
不定積分
問題1
以下の不定積分を求めよ.
$$\int \sin^{-1}xdx$$
解法1
MATLABでの積分はint関数を
答え = int(数式,積分する変数)
と使う. これにより, 以下のように問題を解ける.
syms x % Symbolic変数の宣言
f = asin(x) % 数式の決定
F = int(f,x) % 不定積分の求解
% Visualize
figure
fplot(f)
hold on
fplot(F)
hold off
grid on
legend('Original','Itegration')
saveas(gcf,'int1.png')
f =
asin(x)
F =
x*asin(x) + (1 - x^2)^(1/2)
定積分
問題2
以下の積分を求めよ.
$$\int^1_0\frac{x}{x^4+1}dx$$
解法2
定積分の場合は,
答え = int(数式,積分する変数, 積分区間始点,積分区間終点)
となります. よって,
syms x % Symbolic変数の宣言
f = x/(x^4+1) % 数式の決定
F = int(f,x,0,1) % 定積分
F_for_graph = int(f,x) % グラフ作成用不定積分
% Visualize
figure
fplot(f)
hold on
fplot(F_for_graph)
hold off
grid on
legend('Original','Itegration')
saveas(gcf,'int2.png')
f =
x/(x^4 + 1)
F =
pi/8
F_for_graph =
atan(x^2)/2
重積分
問題3
以下の重積分を計算せよ.
$$ \int^{2\pi}_0 \int^{2\pi}_0 \sin(2x+y) dxdy $$
解法3
syms x y % Symbolic変数の宣言
f = sin(2*x+y) % 数式の決定
F1 = int(f,x,0,pi/2) % xで定積分
F2 = int(F1,y,0,pi/2) % yで定積分
F1_for_graph = int(f,x) % グラフ作成用不定積分1
F2_for_graph = int(F1_for_graph,y) % グラフ作成用不定積分2
% Visualize
figure
subplot(3,1,1)
fsurf(f)
grid on
title('f')
subplot(3,1,2)
fsurf(F1_for_graph)
grid on
title('F_x')
subplot(3,1,3)
fsurf(F2_for_graph)
grid on
title('F_{xy}')
f =
sin(2*x + y)
F1 =
cos(y)
F2 =
1
F1_for_graph =
-cos(2*x + y)/2
F2_for_graph =
-sin(2*x + y)/2
終わりに
今回は微分積分学の基本的な計算をMATLABで行う方法を記載させていただきました.
もし, この記事を見て, 自分の計算を確認する術を知った人が, テストで高得点を取ってくれたらうれしいです!!
参考文献
宮島静雄著, 微分積分学Ⅰ, 共立出版株式会社, 2015
宮島静雄著, 微分積分学Ⅱ, 共立出版株式会社, 2015
三宅敏恒著, 入門微分積分, 培風館, 2016