静岡理工科大学 菅沼ホーム 目次 索引

微分方程式(ルンゲ・クッタ)

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

  プログラムは,以下の微分方程式をルンゲ・クッタ法によって,0 から 1 秒まで解いた例です.

  d2y/dt2 + 3dy/dt + 2y = 1  初期条件はすべて0
  (x[0] = y, x[1] = dy/dt)

  1. C++

    /**********************************/
    /* 微分方程式(ルンゲ・クッタ法) */
    /*      例:d2x/dt+3dx/dt+2x=1    */
    /*      coded by Y.Suganuma       */
    /**********************************/
    #include <stdio.h>
    
    void snx(double, double *, double *);
    double rkg(double, double, double *, double *, double **, int,
               void (*)(double, double *, double *));
    int main()
    {
    	double time, h, *x, *dx, **g;
    	int i1, n = 2;
    /*
              初期設定
    */
    	x  = new double [n];
    	dx = new double [n];
    	g  = new double * [4];
    	for (i1 = 0; i1 < 4; i1++)
    		g[i1] = new double [n];
    	time = 0.0;
    	h    = 0.01;
    	x[0] = 0.0;
    	x[1] = 0.0;
    /*
              計算と出力
    */
    	for (i1 = 0; i1 < 101; i1++) {
    		printf("time = %f, x = %f\n", time, x[0]);
    		time = rkg(time, h, x, dx, g, n, snx);
    	}
    
    	return 0;
    }
    
    /****************/
    /* 微係数の計算 */
    /****************/
    void snx(double time, double *x, double *dx)
    {
    	dx[0] = x[1];
    	dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0;
    }
    
    /*******************************************/
    /* ルンゲ・クッタ法  dx/dt=f(t,x)          */
    /*      time : 現在の時間                  */
    /*      h : 時間刻み幅                     */
    /*      x : 現在の状態                     */
    /*      dx : 微係数(f(t,x):snxで計算する)*/
    /*      g : 作業域(g[4][n])              */
    /*      n : 微分方程式の次数               */
    /*      sub : 微係数を計算する関数の名前   */
    /*      return : time+h                    */
    /*******************************************/
    double rkg(double time, double h, double *x, double *dx, double **g,
               int n, void (*sub)(double, double *, double *))
    {
    	int i1;
    	double h2;
    
    	h2 = 0.5 * h;
    
    	(*sub)(time, x, dx);
    	for (i1 = 0; i1 < n; i1++)
    		g[0][i1] = h * dx[i1];
    
    	time += h2;
    	for (i1 = 0; i1 < n; i1++)
    		g[1][i1] = x[i1] + 0.5 * g[0][i1];
    	(*sub)(time, g[1], dx);
    	for (i1 = 0; i1 < n; i1++)
    		g[1][i1] = h * dx[i1];
    
    	for (i1 = 0; i1 < n; i1++)
    		g[2][i1] = x[i1] + 0.5 * g[1][i1];
    	(*sub)(time, g[2], dx);
    	for (i1 = 0; i1 < n; i1++)
    		g[2][i1] = h * dx[i1];
    
    	time += h2;
    	for (i1 = 0; i1 < n; i1++)
    		g[3][i1] = x[i1] + g[2][i1];
    	(*sub)(time, g[3], dx);
    	for (i1 = 0; i1 < n; i1++)
    		g[3][i1] = h * dx[i1];
    
    	for (i1 = 0; i1 < n; i1++)
    		x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;
    
    	return time;
    }
    			

  2. Java

    /**********************************/
    /* 微分方程式(ルンゲ・クッタ法) */
    /*      例:d2x/dt+3dx/dt+2x=1    */
    /*      coded by Y.Suganuma       */
    /**********************************/
    import java.io.*;
    
    public class Test {
    	public static void main(String args[]) throws IOException
    	{
    		double time, h, x[], dx[], g[][];
    		int i1, n = 2;
    	/*
    	          初期設定
    	*/
    		x    = new double [n];
    		dx   = new double [n];
    		g    = new double [4][n];
    
    		time = 0.0;
    		h    = 0.01;
    		x[0] = 0.0;
    		x[1] = 0.0;
    	/*
    	          計算と出力
    	*/
    		Kansu kn = new Kansu();
    
    		for (i1 = 0; i1 < 101; i1++) {
    			System.out.println("time = " + time + " x = " + x[0]);
    			time = kn.rkg(time, h, x, dx, g, n);
    		}
    	}
    }
    
    /****************/
    /* 微係数の計算 */
    /****************/
    class Kansu extends RKG
    {
    	void snx(double time, double x[], double dx[])
    	{
    		dx[0] = x[1];
    		dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0;
    	}
    }
    
    abstract class RKG
    {
    	/*************************************************/
    	/* ルンゲ・クッタ法  dx/dt=f(t,x)                */
    	/*      time : 現在の時間                        */
    	/*      h : 時間刻み幅                           */
    	/*      x : 現在の状態                           */
    	/*      dx : 微係数(f(t,x):snxで計算する)      */
    	/*      g : 作業域(g[4][n])                    */
    	/*      n : 微分方程式の次数                     */
    	/*      return : time+h                          */
    	/*************************************************/
    
    	abstract void snx(double time, double x[], double dx[]);
    
    	double rkg(double time, double h, double x[], double dx[], double g[][], int n)
    	{
    		int i1;
    		double h2;
    
    		h2 = 0.5 * h;
    
    		snx(time, x, dx);
    		for (i1 = 0; i1 < n; i1++)
    			g[0][i1] = h * dx[i1];
    
    		time += h2;
    		for (i1 = 0; i1 < n; i1++)
    			g[1][i1] = x[i1] + 0.5 * g[0][i1];
    		snx(time, g[1], dx);
    		for (i1 = 0; i1 < n; i1++)
    			g[1][i1] = h * dx[i1];
    
    		for (i1 = 0; i1 < n; i1++)
    			g[2][i1] = x[i1] + 0.5 * g[1][i1];
    		snx(time, g[2], dx);
    		for (i1 = 0; i1 < n; i1++)
    			g[2][i1] = h * dx[i1];
    
    		time += h2;
    		for (i1 = 0; i1 < n; i1++)
    			g[3][i1] = x[i1] + g[2][i1];
    		snx(time, g[3], dx);
    		for (i1 = 0; i1 < n; i1++)
    			g[3][i1] = h * dx[i1];
    
    		for (i1 = 0; i1 < n; i1++)
    			x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;
    
    		return time;
    	}
    }
    			

  3. JavaScript

      ここをクリックすると,JavaScript の仕様に適合した形で微分方程式を入力することによって,任意の微分方程式の解を画面上で求めることができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>微分方程式(ルンゲ・クッタ)</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		eq = new Array();
    
    		function main()
    		{
    					// データの設定
    			let n  = parseInt(document.getElementById("order").value);
    			let t1 = parseFloat(document.getElementById("t1").value);
    			let t2 = parseFloat(document.getElementById("t2").value);
    			let h  = parseFloat(document.getElementById("haba").value);
    			let x  = new Array();
    			x = document.getElementById("ini").value.split(/ {1,}/)
    			for (let i1 = 0; i1 < n; i1++)
    				x[i1] = parseFloat(x[i1]);
    			eq = document.getElementById("func").value.split(/\n{1,}/)
    			for (let i1 = 0; i1 < n; i1++)
    				eq[i1] = eq[i1] + ";";
    			dx = new Array();
    			g  = new Array();
    			for (let i1 = 0; i1 < 4; i1++)
    				g[i1] = new Array();
    					// 実行,結果
    			let time = t1;
    			let c = 100000;
    			let str  = "";
    			while (time < t2+0.1*h) {
    				let s1 = Math.round(time * c) / c;
    				str = str + "time = " + s1 + " x =";
    				for (let i1 = 0; i1 < n; i1++) {
    					let s2 = Math.round(x[i1] * c) / c;
    					str = str + " " + s2;
    				}
    				str += "\n";
    				time = rkg(time, h, x, dx, g, n, snx);
    			}
    			document.getElementById("res").value = str;
    		}
    
    		/****************/
    		/* 関数値の計算 */
    		/****************/
    		function snx(n, x, dx)
    		{
    			for (let i1 = 0; i1 < n; i1++)
    				dx[i1] = eval(eq[i1]);
    		}
    
    		/********************************************/
    		/* ルンゲ・クッタ法  dx/dt=f(t,x)           */
    		/*      time : 現在の時間                   */
    		/*      h : 時間刻み幅                      */
    		/*      x : 現在の状態                      */
    		/*      dx : 微係数(f(t,x):snxで計算する) */
    		/*      g : 作業域(g[4][n])               */
    		/*      n : 微分方程式の次数                */
    		/*      fn : 関数値を計算する関数      */
    		/*      return : time+h                     */
    		/********************************************/
    		function rkg(time, h, x, dx, g, n, fn)
    		{
    			let h2 = 0.5 * h;
    
    			fn(n, x, dx);
    			for (let i1 = 0; i1 < n; i1++)
    				g[0][i1] = h * dx[i1];
    
    			time += h2;
    			for (let i1 = 0; i1 < n; i1++)
    				g[1][i1] = x[i1] + 0.5 * g[0][i1];
    			fn(n, g[1], dx);
    			for (let i1 = 0; i1 < n; i1++)
    				g[1][i1] = h * dx[i1];
    
    			for (let i1 = 0; i1 < n; i1++)
    				g[2][i1] = x[i1] + 0.5 * g[1][i1];
    			fn(n, g[2], dx);
    			for (let i1 = 0; i1 < n; i1++)
    				g[2][i1] = h * dx[i1];
    
    			time += h2;
    			for (let i1 = 0; i1 < n; i1++)
    				g[3][i1] = x[i1] + g[2][i1];
    			fn(n, g[3], dx);
    			for (let i1 = 0; i1 < n; i1++)
    				g[3][i1] = h * dx[i1];
    
    			for (let i1 = 0; i1 < n; i1++)
    				x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;
    
    			return time;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>微分方程式(ルンゲ・クッタ)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールド及びテキストエリアには,例として,微分方程式,
    		<P STYLE="text-align:center"><IMG SRC="rkg.gif"></P>
    		<DT>を解く場合に対する値が設定されています(他の問題を実行する場合は,それらを適切に修正してください).ルンゲ・クッタ法を使用するためには,1 次の連立微分方程式に変換する必要があるため,
    		<P STYLE="text-align:center"><IMG SRC="rkg1.gif"></P>
    		<DT>とおき,以下に示すような 1 次の連立微分方程式に変換しています.
    		<P STYLE="text-align:center"><IMG SRC="rkg2.gif"></P>
    		<DT>テキストエリアに設定してあるのは,これらの式の右辺です.ただし,変数 x<SUB>1</SUB> や x<SUB>2</SUB> の代わりに,配列の要素 x[0] と x[1] を使用しています.なお,式は,JavaScript の仕様に適合した形式で記述してあることに注意してください.
    	</DL>
    
    	<DIV STYLE="text-align:center">
    
    		次数(n):<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="2"> 
    		範囲:<INPUT ID="t1" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.0">,
    				<INPUT ID="t2" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="1.0"> 
    		幅:<INPUT ID="haba" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.01"> 
    		初期値(n個):<INPUT ID="ini" STYLE="font-size: 100%" TYPE="text" SIZE="10" VALUE="0.0 0.0"><BR><BR>
    		微分方程式  <BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR>
    		dx = <TEXTAREA ID="func" STYLE="font-size: 110%" COLS="70" ROWS="10">
    x[1]
    -2.0 * x[0] - 3.0 * x[1] + 1.0
    		</TEXTAREA><BR><BR>
    		結果<BR>
    		<TEXTAREA ID="res" STYLE="font-size: 110%" COLS="50" ROWS="10"></TEXTAREA>
    	</DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /*
              初期設定
    */
    	$n  = 2;
    	$x  = array($n);
    	$dx = array($n);
    	$g  = array(4);
    	for ($i1 = 0; $i1 < 4; $i1++)
    		$g[$i1] = array($n);
    	$time = 0.0;
    	$h    = 0.01;
    	$x[0] = 0.0;
    	$x[1] = 0.0;
    /*
              計算と出力
    */
    	for ($i1 = 0; $i1 < 101; $i1++) {
    		printf("time = %f, x = %f\n", $time, $x[0]);
    		$time = rkg($time, $h, $x, $dx, $g, $n, "snx");
    	}
    
    /****************/
    /* 微係数の計算 */
    /****************/
    function snx($time, $x, &$dx)
    {
    	$dx[0] = $x[1];
    	$dx[1] = -2.0 * $x[0] - 3.0 * $x[1] + 1.0;
    }
    
    /*******************************************/
    /* ルンゲ・クッタ法  dx/dt=f(t,x)          */
    /*      time : 現在の時間                  */
    /*      h : 時間刻み幅                     */
    /*      x : 現在の状態                     */
    /*      dx : 微係数(f(t,x):snxで計算する)*/
    /*      g : 作業域(g[4][n])              */
    /*      n : 微分方程式の次数               */
    /*      sub : 微係数を計算する関数の名前   */
    /*      return : time+h                    */
    /*******************************************/
    function rkg($time, $h, &$x, $dx, $g, $n, $sub)
    {
    	$h2 = 0.5 * $h;
    
    	$sub($time, $x, $dx);
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$g[0][$i1] = $h * $dx[$i1];
    
    	$time += $h2;
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$g[1][$i1] = $x[$i1] + 0.5 * $g[0][$i1];
    	$sub($time, $g[1], $dx);
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$g[1][$i1] = $h * $dx[$i1];
    
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$g[2][$i1] = $x[$i1] + 0.5 * $g[1][$i1];
    	$sub($time, $g[2], $dx);
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$g[2][$i1] = $h * $dx[$i1];
    
    	$time += $h2;
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$g[3][$i1] = $x[$i1] + $g[2][$i1];
    	$sub($time, $g[3], $dx);
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$g[3][$i1] = $h * $dx[$i1];
    
    	for ($i1 = 0; $i1 < $n; $i1++)
    		$x[$i1] = $x[$i1] + ($g[0][$i1] + 2.0 * $g[1][$i1] + 2.0 * $g[2][$i1] + $g[3][$i1]) / 6.0;
    
    	return $time;
    }
    
    ?>
    			

  5. Ruby

    #*********************************/
    # 微分方程式(ルンゲ・クッタ法) */
    #      例:d2x/dt+3dx/dt+2x=1    */
    #      coded by Y.Suganuma       */
    #*********************************/
    
    #***************/
    # 微係数の計算 */
    #***************/
    snx = Proc.new { |time, x, dx|
    	dx[0] = x[1]
    	dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0
    }
    
    #******************************************/
    # ルンゲ・クッタ法  dx/dt=f(t,x)          */
    #      time : 現在の時間                  */
    #      h : 時間刻み幅                     */
    #      x : 現在の状態                     */
    #      dx : 微係数(f(t,x):snxで計算する)*/
    #      g : 作業域(g[4][n])              */
    #      n : 微分方程式の次数               */
    #      snx : 微係数を計算する関数の名前   */
    #      return : time+h                    */
    #******************************************/
    def rkg(time, h, x, dx, g, n, &sub)
    
    	h2 = 0.5 * h
    
    	sub.call(time, x, dx)
    	for i1 in 0 ... n
    		g[0][i1] = h * dx[i1]
    	end
    
    	time += h2
    	for i1 in 0 ... n
    		g[1][i1] = x[i1] + 0.5 * g[0][i1]
    	end
    	sub.call(time, g[1], dx)
    	for i1 in 0 ... n
    		g[1][i1] = h * dx[i1]
    	end
    
    	for i1 in 0 ... n
    		g[2][i1] = x[i1] + 0.5 * g[1][i1]
    	end
    	sub.call(time, g[2], dx)
    	for i1 in 0 ... n
    		g[2][i1] = h * dx[i1]
    	end
    
    	time += h2
    	for i1 in 0 ... n
    		g[3][i1] = x[i1] + g[2][i1]
    	end
    	sub.call(time, g[3], dx)
    	for i1 in 0 ... n
    		g[3][i1] = h * dx[i1]
    	end
    
    	for i1 in 0 ... n
    		x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0
    	end
    
    	return time
    end
    
    		# 初期設定
    n  = 2
    x  = Array.new(n)
    dx = Array.new(n)
    g  = Array.new(4)
    for i1 in 0 ... 4
    	g[i1] = Array.new(n)
    end
    time = 0.0
    h    = 0.01
    x[0] = 0.0
    x[1] = 0.0
    		# 計算と出力
    for i1 in 0 ... 101
    	printf("time = %f, x = %f\n", time, x[0])
    	time = rkg(time, h, x, dx, g, n, &snx)
    end
    			

  6. Python

    # -*- coding: UTF-8 -*-
    import numpy as np
    from math import *
    
    ############################################
    # ルンゲ・クッタ法  dx/dt=f(t,x)
    #      time : 現在の時間
    #      h : 時間刻み幅
    #      x : 現在の状態
    #      dx : 微係数(f(t,x):subで計算する)
    #      g : 作業域(g[4][n])
    #      n : 微分方程式の次数
    #      sub : 微係数を計算する関数の名前
    #      return : time+h
    #      coded by Y.Suganuma
    ############################################
    
    def rkg(time, h, x, dx, g, n, sub) :
    
    	h2 = 0.5 * h
    	sub(time, x, dx)
    	for i1 in range(0, n) :
    		g[0][i1] = h * dx[i1]
    
    	time += h2
    	for i1 in range(0, n) :
    		g[1][i1] = x[i1] + 0.5 * g[0][i1]
    	sub(time, g[1], dx)
    	for i1 in range(0, n) :
    		g[1][i1] = h * dx[i1]
    
    	for i1 in range(0, n) :
    		g[2][i1] = x[i1] + 0.5 * g[1][i1]
    	sub(time, g[2], dx)
    	for i1 in range(0, n) :
    		g[2][i1] = h * dx[i1]
    
    	time += h2
    	for i1 in range(0, n) :
    		g[3][i1] = x[i1] + g[2][i1]
    	sub(time, g[3], dx)
    	for i1 in range(0, n) :
    		g[3][i1] = h * dx[i1]
    
    	for i1 in range(0, n) :
    		x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0
    
    	return time
    
    ############################################
    # 微分方程式(ルンゲ・クッタ法)
    #      例:d2x/dt+3dx/dt+2x=1
    #      coded by Y.Suganuma
    ############################################
    
    			# 微係数の計算 */
    def snx(time, x, dx) :
    	dx[0] = x[1]
    	dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0
    			# 初期設定
    n    = 2
    x    = np.array([0.0, 0.0], np.float)
    dx   = np.empty(n, np.float)
    g    = np.empty((4, n), np.float)
    time = 0.0
    h    = 0.01
    			# 計算と出力
    for i1 in range(0, 101) :
    	print("time = " + str(time) + ", x = " + str(x[0]))
    	time = rkg(time, h, x, dx, g, n, snx)
    			

  7. C#

    /**********************************/
    /* 微分方程式(ルンゲ・クッタ法) */
    /*      例:d2x/dt+3dx/dt+2x=1    */
    /*      coded by Y.Suganuma       */
    /**********************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	public Test1()
    	{
    	/*
    	          初期設定
    	*/
    		int n        = 2;
    		double time  = 0.0;
    		double h     = 0.01;
    		double[] x   = {0.0, 0.0};
    		double[] dx  = new double [n];
    		double[][] g = new double [4][];
    		g[0] = new double[2];
    		g[1] = new double[2];
    		g[2] = new double[2];
    		g[3] = new double[2];
    	/*
    	          計算と出力
    	*/
    		for (int i1 = 0; i1 < 101; i1++) {
    			Console.WriteLine("time = " + time + " x = " + x[0]);
    			time = rkg(time, h, x, dx, g, n, snx);
    		}
    	}
    			// 微係数の計算
    	int snx(double time, double[] x, double[] dx)
    	{
    		dx[0] = x[1];
    		dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0;
    		return 0;
    	}
    
    	/********************************************/
    	/* ルンゲ・クッタ法  dx/dt=f(t,x)           */
    	/*      time : 現在の時間                   */
    	/*      h : 時間刻み幅                      */
    	/*      x : 現在の状態                      */
    	/*      dx : 微係数(f(t,x):snxで計算する) */
    	/*      g : 作業域(g[4][n])               */
    	/*      n : 微分方程式の次数                */
    	/*      fn : 関数値を計算する関数           */
    	/*      return : time+h                     */
    	/********************************************/
    	double rkg(double time, double h, double[] x, double[] dx, double[][] g,
    	           int n, Func<double, double[], double[], int> fn)
    	{
    		double h2 = 0.5 * h;
    
    		fn(time, x, dx);
    		for (int i1 = 0; i1 < n; i1++)
    			g[0][i1] = h * dx[i1];
    
    		time += h2;
    		for (int i1 = 0; i1 < n; i1++)
    			g[1][i1] = x[i1] + 0.5 * g[0][i1];
    		fn(time, g[1], dx);
    		for (int i1 = 0; i1 < n; i1++)
    			g[1][i1] = h * dx[i1];
    
    		for (int i1 = 0; i1 < n; i1++)
    			g[2][i1] = x[i1] + 0.5 * g[1][i1];
    		fn(time, g[2], dx);
    		for (int i1 = 0; i1 < n; i1++)
    			g[2][i1] = h * dx[i1];
    
    		time += h2;
    		for (int i1 = 0; i1 < n; i1++)
    			g[3][i1] = x[i1] + g[2][i1];
    		fn(time, g[3], dx);
    		for (int i1 = 0; i1 < n; i1++)
    			g[3][i1] = h * dx[i1];
    
    		for (int i1 = 0; i1 < n; i1++)
    			x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;
    
    		return time;
    	}
    }
    			

  8. VB

    ''''''''''''''''''''''''''''''''''
    ' 微分方程式(ルンゲ・クッタ法) '
    '      例:d2x/dt+3dx/dt+2x=1    '
    '      coded by Y.Suganuma       '
    ''''''''''''''''''''''''''''''''''
    Module Test
    
    	Sub Main()
    				'
    				' 初期設定
    				'
    		Dim n As Integer   = 2
    		Dim time As Double = 0.0
    		Dim h As Double    = 0.01
    		Dim x() As Double  = {0.0, 0.0}
    		Dim dx(n) As Double
    		Dim g0(2) As Double
    		Dim g1(2) As Double
    		Dim g2(2) As Double
    		Dim g3(2) As Double
    				' 微係数の計算
    		Dim snx = Function(tm, v, dv)
    			dv(0) = v(1)
    			dv(1) = -2.0 * v(0) - 3.0 * v(1) + 1.0
    			Return 0
    		End Function
    				'
    				' 計算と出力
    				'
    		For i1 As Integer = 0 To 100
    			Console.WriteLine("time = " & time & " x = " & x(0))
    			time = rkg(time, h, x, dx, g0, g1, g2, g3, n, snx)
    		Next
    
    	End Sub
    
    	''''''''''''''''''''''''''''''''''''''''''''
    	' ルンゲ・クッタ法  dx/dt=f(t,x)           '
    	'      time : 現在の時間                   '
    	'      h : 時間刻み幅                      '
    	'      x : 現在の状態                      '
    	'      dx : 微係数(f(t,x):snxで計算する) '
    	'      g0,g1,g2,g3 : 作業域(gi(n))       '
    	'      n : 微分方程式の次数                '
    	'      fn : 関数値を計算する関数           '
    	'      return : time+h                     '
    	''''''''''''''''''''''''''''''''''''''''''''
    	Function rkg(time As Double, h As Double, x() As Double, dx() As Double,
    	             g0() As Double, g1() As Double, g2() As Double, g3() As Double,
    	             n As Integer, fn As Func(Of Double, Double(), Double(), Integer))
    
    		Dim h2 As Double = 0.5 * h
    
    		fn(time, x, dx)
    		For i1 As Integer = 0 To n-1
    			g0(i1) = h * dx(i1)
    		Next
    
    		time += h2
    		For i1 As Integer = 0 To n-1
    			g1(i1) = x(i1) + 0.5 * g0(i1)
    		Next
    		fn(time, g1, dx)
    		For i1 As Integer = 0 To n-1
    			g1(i1) = h * dx(i1)
    		Next
    
    		For i1 As Integer = 0 To n-1
    			g2(i1) = x(i1) + 0.5 * g1(i1)
    		Next
    		fn(time, g2, dx)
    		For i1 As Integer = 0 To n-1
    			g2(i1) = h * dx(i1)
    		Next
    
    		time += h2
    		For i1 As Integer = 0 To n-1
    			g3(i1) = x(i1) + g2(i1)
    		Next
    		fn(time, g3, dx)
    		For i1 As Integer = 0 To n-1
    			g3(i1) = h * dx(i1)
    		Next
    
    		For i1 As Integer = 0 To n-1
    			x(i1) = x(i1) + (g0(i1) + 2.0 * g1(i1) + 2.0 * g2(i1) + g3(i1)) / 6.0
    		Next
    
    		Return time
    
    	End Function
    
    End Module
    			

静岡理工科大学 菅沼ホーム 目次 索引