演習
演習 1
のチャートを、の範囲で幅 300、高さ 300 のウィンドウに描画しましょう。
演習 2
についてのチャートを、、の範囲で幅 300、高さ 300 のウィンドウに描画しましょう。
演習 3
以下の初期値問題について、のチャートを、の範囲で幅 300、高さ 300 のウィンドウに描画しましょう。
演習 4
ある時点 において、細胞の突然変異率 の変化を照射される放射線の線量率 によって以下の式で表すとします。
は以下の値の定数です。
| 記号 | 値 |
|---|---|
初期値 を次式としたとき、 のチャートを の範囲で描画しましょう。
なお、の範囲はに対して非常に小さいので、lineの呼び出し時に縦軸の値を 倍すると良い。
演習 5
次の初期値問題の解析解を求めましょう。
そして、解析解と Euler 法による数値解をそれぞれ、、の範囲で幅 300、高さ 300 のウィンドウに描画しましょう。
演習 6
次の初期値問題の Euler 法による数値解を 、の範囲で幅 300、高さ 300 のウィンドウに描画しましょう。
なお、この連立微分方程式は Lotka-Volterra の方程式 として知られています。 Lotka-Volterra の方程式は生物の捕食被食関係による個体数の増減を表現する数理モデルです。
演習 7
演習 5 を Runge-Kutta 法を用いて解き、Euler 法での結果と比較しましょう。
演習 8
演習 6 を Runge-Kutta 法を用いて解き、Euler 法での結果と比較しましょう。
演習 9
2 次の Runge-Kutta 法（Runge-Kutta の中点法、Heun 法、改良オイラー法）を用いて次の初期値問題を解きましょう。
演習 10
以下の連立微分方程式の初期値問題を考えます。
から始まり、球の中心の x 座標を 、y 座標を とした時の球の運動のアニメーションを、、の範囲で幅 300、高さ 300 のウィンドウに描画しましょう。
演習 11
2 つの天体が座標 と に固定されていて、それぞれの質量を と とします。 これらの天体の周りを巡る質量 の衛星の運動方程式は以下の微分方程式で表されます。
、 とし、新たな未知関数 と を導入すると、以下のように書き直すことができます。
以下の初期値が与えられたときの衛星の運動のアニメーションを、、の範囲で幅 300、高さ 300 のウィンドウに描画しましょう。
演習 12
二次元空間上のある点 の速度 が次式で与えられている流れ場について、 を自由に決め、流線を描いてみましょう。
演習 13
以下のプログラムの vx と vy の値を自由に決め、流線を描いてみましょう。
const xLeft = 0;
const xRight = 4;
const yTop = 4;
const yBottom = 0;
let t1 = 0;
const x1 = [2.5, 2.5];
function setup() {
createCanvas(300, 300);
background(255);
applyTransform();
drawVectorField(20);
}
function draw() {
if (x1[0] < xLeft || xRight < x1[0] || x1[1] < yBottom || yTop < x1[1]) {
return;
}
applyTransform();
const h = 0.1;
const d = x1.length;
const x2 = new Array(d);
const tmp = new Array(d);
const f1 = new Array(d);
const f2 = new Array(d);
const f3 = new Array(d);
const f4 = new Array(d);
const t2 = t1 + h;
dx(t1, x1, f1);
for (let j = 0; j < d; ++j) {
tmp[j] = x1[j] + (h * f1[j]) / 2;
}
dx(t1, tmp, f2);
for (let j = 0; j < d; ++j) {
tmp[j] = x1[j] + (h * f2[j]) / 2;
}
dx(t1, tmp, f3);
for (let j = 0; j < d; ++j) {
tmp[j] = x1[j] + h * f3[j];
}
dx(t1, tmp, f4);
for (let j = 0; j < d; ++j) {
x2[j] = x1[j] + (h * (f1[j] + 2 * f2[j] + 2 * f3[j] + f4[j])) / 6;
}
stroke(0, 0, 255);
line(x1[0], x1[1], x2[0], x2[1]);
t1 = t2;
for (let j = 0; j < d; ++j) {
x1[j] = x2[j];
}
}
function dx(t, xt, dxt) {
const vx = [
[-2.0, -1.0, 0.0, 1.0, 2.0],
[-2.0, -1.0, 0.0, 1.0, 2.0],
[-2.0, -1.0, 0.0, 1.0, 2.0],
[-2.0, -1.0, 0.0, 1.0, 2.0],
[-2.0, -1.0, 0.0, 1.0, 2.0],
];
const vy = [
[2.0, 2.0, 2.0, 2.0, 2.0],
[1.0, 1.0, 1.0, 1.0, 1.0],
[0.0, 0.0, 0.0, 0.0, 0.0],
[-1.0, -1.0, -1.0, -1.0, -1.0],
[-2.0, -2.0, -2.0, -2.0, -2.0],
];
const i = constrain(int(xt[1]), 0, 3);
const j = constrain(int(xt[0]), 0, 3);
const vx1 = (xt[0] - i) * vx[i + 1][j] - (xt[0] - (i + 1)) * vx[i][j];
const vx2 = (xt[0] - i) * vx[i + 1][j + 1] - (xt[0] - (i + 1)) * vx[i][j + 1];
dxt[0] = (xt[1] - j) * vx2 - (xt[1] - (j + 1)) * vx1;
const vy1 = (xt[0] - i) * vy[i + 1][j] - (xt[0] - (i + 1)) * vy[i][j];
const vy2 = (xt[0] - i) * vy[i + 1][j + 1] - (xt[0] - (i + 1)) * vy[i][j + 1];
dxt[1] = (xt[1] - j) * vy2 - (xt[1] - (j + 1)) * vy1;
}
function applyTransform() {
resetMatrix();
scale(width / (xRight - xLeft), height / (yBottom - yTop));
translate(-xLeft, -yTop);
strokeWeight((xRight - xLeft) / width);
}
function drawVectorField(steps) {
const tmpX = new Array(2);
const tmpDx = new Array(2);
const xStep = (xRight - xLeft) / steps;
const yStep = (yTop - yBottom) / steps;
for (let x = xLeft; x <= xRight; x += xStep) {
for (let y = yBottom; y <= yTop; y += yStep) {
tmpX[0] = x;
tmpX[1] = y;
dx(0, tmpX, tmpDx);
drawVector(tmpX[0], tmpX[1], tmpDx[0], tmpDx[1]);
}
}
}
function drawVector(x, y, dx, dy) {
const d = dist(0, 0, dx, dy) / 50;
const theta = atan2(dy, dx);
stroke(0);
line(x, y, d * cos(theta) + x, d * sin(theta) + y);
}
function mouseClicked() {
x1[0] = ((xRight - xLeft) * float(mouseX)) / width + xLeft;
x1[1] = ((yBottom - yTop) * float(mouseY)) / height + yTop;
}