{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2022年度プログラミング演習A・B\n",
"\n",
"# 第14回レポート課題の解説"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習1\n",
"\n",
"連立一次方程式を解くことによって、関数\n",
"\n",
"$$\n",
"f(x)=e^x\n",
"$$\n",
"\n",
"の節点\n",
"\n",
"$$\n",
"x_0=0,\\quad x_1=1\n",
"$$\n",
"\n",
"におけるエルミート補間多項式を求め、$-2\\leq x\\leq 3$ の範囲で元の関数 $f(x)$ とエルミート補間多項式のグラフを描画してください。\n",
"\n",
"なお、連立一次方程式を解く際には演算子`\\`を使用し、\n",
"\n",
"$$\n",
"f'(x)=e^x\n",
"$$\n",
"\n",
"は既知のこととしてよいです。"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"\n",
"function value = f(x) %関数f(x)\n",
" value = exp(x);\n",
"end\n",
"\n",
"function value = f_dif(x) %関数f(x)の微分\n",
" value = exp(x);\n",
"end\n",
"\n",
"x = -2:0.01:3; %x座標の配列\n",
"y = f(x);\n",
"plot(x,y,\"DisplayName\",\"f(x)\") %f(x)のグラフを描画\n",
"grid on\n",
"hold on\n",
"\n",
"x0 = 0; %節点\n",
"x1 = 1;\n",
"\n",
"A = [1,x0,x0^2,x0^3; 1,x1,x1^2,x1^3; 0,1,2*x0,3*x0^2; 0,1,2*x1,3*x1^2];\n",
"b = [f(x0); f(x1); f_dif(x0); f_dif(x1)];\n",
"c = A\\b; %連立一次方程式Ac=bを解いて係数ベクトルcを求める\n",
"\n",
"y = c(1)+c(2)*x+c(3)*x.^2+c(4)*x.^3; %エルミート補間多項式を計算\n",
"plot(x,y,\"DisplayName\",\"Hermite interpolating polynomial\") %エルミート補間多項式のグラフを描画\n",
"\n",
"legend %凡例の表示"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### （別解）"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"\n",
"function value = f(x) %関数f(x)\n",
" value = exp(x);\n",
"end\n",
"\n",
"function value = f_dif(x) %関数f(x)の微分\n",
" value = exp(x);\n",
"end\n",
"\n",
"x = -2:0.01:3; %x座標の配列\n",
"y = f(x);\n",
"plot(x,y,\"DisplayName\",\"f(x)\") %f(x)のグラフを描画\n",
"grid on\n",
"hold on\n",
"\n",
"node = [0;1]; %節点の配列（縦ベクトル）\n",
"n = length(node)-1; %nodeの要素数が変わっても対応できるように、以下ではnを使う\n",
"\n",
"A = zeros(2*n+2);\n",
"for j = 0:2*n+1\n",
" A(1:n+1,j+1) = node.^j; %行列Aの上半分を列ごとに計算\n",
"end\n",
"for j = 1:2*n+1\n",
" A(n+2:2*n+2,j+1) = j*node.^(j-1); %行列Aの下半分を列ごとに計算\n",
"end\n",
"b = [f(node);f_dif(node)];\n",
"c = A\\b; %連立一次方程式Ac=bを解いて係数ベクトルcを求める\n",
"\n",
"y = 0;\n",
"for j = 0:2*n+1\n",
" y = y+c(j+1)*x.^j; %エルミート補間多項式を計算\n",
"end\n",
"plot(x,y,\"DisplayName\",\"Hermite interpolating polynomial\") %エルミート補間多項式のグラフを描画\n",
"\n",
"legend %凡例の表示"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習2\n",
"\n",
"エルミート基底多項式を用いることによって、関数\n",
"\n",
"$$\n",
"f(x)=e^x\n",
"$$\n",
"\n",
"の節点\n",
"\n",
"$$\n",
"x_0=0,\\quad x_1=1\n",
"$$\n",
"\n",
"におけるエルミート補間多項式を求め、$-2\\leq x\\leq 3$ の範囲で元の関数 $f(x)$ とエルミート補間多項式のグラフを描画してください。\n",
"\n",
"なお、\n",
"\n",
"$$\n",
"f'(x)=e^x\n",
"$$\n",
"\n",
"は既知のこととしてよいです。"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"\n",
"function value = f(x) %関数f(x)\n",
" value = exp(x);\n",
"end\n",
"\n",
"function value = f_dif(x) %関数f(x)の微分\n",
" value = exp(x);\n",
"end\n",
"\n",
"x = -2:0.01:3; %x座標の配列\n",
"y = f(x);\n",
"plot(x,y,\"DisplayName\",\"f(x)\") %f(x)のグラフを描画\n",
"grid on\n",
"hold on\n",
"\n",
"x0 = 0; %節点\n",
"x1 = 1;\n",
"\n",
"H0 = (1-2*(x-x0)/(x0-x1)).*((x-x1)/(x0-x1)).^2; %エルミート基底多項式を計算\n",
"H1 = (1-2*(x-x1)/(x1-x0)).*((x-x0)/(x1-x0)).^2;\n",
"K0 = (x-x0).*((x-x1)/(x0-x1)).^2;\n",
"K1 = (x-x1).*((x-x0)/(x1-x0)).^2;\n",
"\n",
"y = f(x0)*H0+f(x1)*H1+f_dif(x0)*K0+f_dif(x1)*K1; %エルミート補間多項式を計算\n",
"plot(x,y,\"DisplayName\",\"Hermite interpolating polynomial\") %エルミート補間多項式のグラフを描画\n",
"\n",
"legend %凡例の表示"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### （別解）\n",
"\n",
"ラグランジュ基底多項式\n",
"\n",
"$$\n",
"L_j(x)=\\prod_{k=0\\\\ k\\neq j}^n\\frac{x-x_k}{x_j-x_k}\\quad (j=0,\\ldots,n)\n",
"$$\n",
"\n",
"の微分が\n",
"\n",
"$$\n",
"L_j'(x)=L_j(x)\\sum_{k=0\\\\ k\\neq j}^n\\frac{1}{x-x_k}\\quad (j=0,\\ldots,n)\n",
"$$\n",
"\n",
"と表されること（対数微分により導出可能）を用いれば、エルミート基底多項式\n",
"\n",
"$$\n",
"H_j(x)=\\left(1-2L_j'(x_j)(x-x_j)\\right)\\left(L_j(x)\\right)^2\\quad (j=0,\\ldots,n),\n",
"$$\n",
"\n",
"$$\n",
"K_j(x)=(x-x_j)\\left(L_j(x)\\right)^2\\quad (j=0,\\ldots,n)\n",
"$$\n",
"\n",
"を一般の場合にも計算できる。"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"\n",
"function value = f(x) %関数f(x)\n",
" value = exp(x);\n",
"end\n",
"\n",
"function value = f_dif(x) %関数f(x)の微分\n",
" value = exp(x);\n",
"end\n",
"\n",
"function value = Lagrange(node,j,x) %ラグランジュ基底多項式\n",
" n = length(node)-1;\n",
" numerator = 1;\n",
" denominator = 1;\n",
" for k = [0:j-1,j+1:n]\n",
" numerator = numerator.*(x-node(k+1));\n",
" denominator = denominator*(node(j+1)-node(k+1));\n",
" end\n",
" value = numerator/denominator;\n",
"end\n",
"\n",
"function value = Lagrange_dif(node,j,x) %ラグランジュ基底多項式の微分\n",
" n = length(node)-1;\n",
" tmp = 0;\n",
" for k = [0:j-1,j+1:n]\n",
" tmp = tmp+1./(x-node(k+1));\n",
" end\n",
" value = tmp.*Lagrange(node,j,x);\n",
"end\n",
"\n",
"function value = H(node,j,x) %エルミート基底多項式1\n",
" value = (1-2*Lagrange_dif(node,j,node(j+1))*(x-node(j+1))).*Lagrange(node,j,x).^2;\n",
"end\n",
"\n",
"function value = K(node,j,x) %エルミート基底多項式2\n",
" value = (x-node(j+1)).*Lagrange(node,j,x).^2;\n",
"end\n",
"\n",
"x = -2:0.01:3; %x座標の配列\n",
"y = f(x);\n",
"plot(x,y,\"DisplayName\",\"f(x)\") %f(x)のグラフを描画\n",
"grid on\n",
"hold on\n",
"\n",
"node = [0,1]; %節点の配列\n",
"n = length(node)-1;\n",
"\n",
"y = 0;\n",
"for j = 0:n\n",
" y = y+f(node(j+1))*H(node,j,x); %エルミート補間多項式を計算\n",
" y = y+f_dif(node(j+1))*K(node,j,x);\n",
"end\n",
"plot(x,y,\"DisplayName\",\"Hermite interpolating polynomial\") %エルミート補間多項式のグラフを描画\n",
"\n",
"legend %凡例の表示"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習3（オプション）\n",
"\n",
"$n$ を自然数とし、関数\n",
"\n",
"$$\n",
"f(x)=\\sin\\left(2\\pi x^2\\right)\n",
"$$\n",
"\n",
"と節点\n",
"\n",
"$$\n",
"x_i=\\frac{i}{2n}\\quad (i=0,\\ldots,2n)\n",
"$$\n",
"\n",
"に対して三角関数による補間を考えます。\n",
"\n",
"$n=1,2,3$ の各場合に補間関数を求め、$0\\leq x\\leq 1$ の範囲で元の関数 $f(x)$ と $3$ 個の補間関数のグラフを描画してください。"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"\n",
"function value = f(x) %関数f(x)\n",
" value = sin(2*pi*x.^2);\n",
"end\n",
"\n",
"x = 0:0.01:1; %x座標の配列\n",
"y = f(x);\n",
"plot(x,y,\"DisplayName\",\"f(x)\") %f(x)のグラフを描画\n",
"grid on\n",
"hold on\n",
"\n",
"for n = [1,2,3] %for文でnを変化させる\n",
" node = (0:1/(2*n):1)'; %節点の配列（縦ベクトル）\n",
" \n",
" A = zeros(2*n+1);\n",
" A(:,1) = 1;\n",
" for j = 1:n\n",
" A(:,j+1) = sin(j*pi*node); %行列Aの各列を計算\n",
" A(:,n+j+1) = cos(j*pi*node);\n",
" end\n",
" b = f(node);\n",
" c = A\\b; %連立一次方程式Ac=bを解いて係数ベクトルcを求める\n",
" \n",
" y = c(1);\n",
" for j = 1:n\n",
" y = y+c(j+1)*sin(j*pi*x); %補間関数を計算\n",
" y = y+c(n+j+1)*cos(j*pi*x);\n",
" end\n",
" plot(x,y,\"DisplayName\",[\"n=\",num2str(n)]) %補間関数のグラフを描画\n",
"end\n",
"\n",
"legend %凡例の表示"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Octave",
"language": "octave",
"name": "octave"
},
"language_info": {
"file_extension": ".m",
"help_links": [
{
"text": "GNU Octave",
"url": "https://www.gnu.org/software/octave/support.html"
},
{
"text": "Octave Kernel",
"url": "https://github.com/Calysto/octave_kernel"
},
{
"text": "MetaKernel Magics",
"url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
}
],
"mimetype": "text/x-octave",
"name": "octave",
"version": "5.2.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}