数値計算@常微分方程式@その5@連立常微分方程式

$k_{1}(n)=f(x_n , y_n )$ (1-1)
$k_{2}(n)=f(x_n+\frac{h}{2},y_n+\frac{h}{2}*k_{1}(n))$ (1-2)
$k_{3}(n)=f(x_n+\frac{h}{2},y_n+\frac{h}{2}*k_{2}(n))$ (1-3)
$k_{4}(n)=f(x_n+h,y_n + k_{3}(n)*h)$ (1-4)

$y_{n+1} = y_n + \frac{1}{6}( k_{1}(n) + 2k_{2}(n) + 2k_{3}(n) + k_{4}(n) )*h$ (1-5)

ここで、を外に出すと

$k_{1}(n)=h*f(x_n , y_n )$ (1-1)
$k_{2}(n)=h*f(x_n+\frac{h}{2},y_n+\frac{1}{2}*k_{1}(n))$ (1-2)
$k_{3}(n)=h*f(x_n+\frac{h}{2},y_n+\frac{1}{2}*k_{2}(n))$ (1-3)
$k_{4}(n)=h*f(x_n+h,y_n + k_{3}(n))$ (1-4)

$y_{n+1} = y_n + \frac{1}{6}( k_{1}(n) + 2k_{2}(n) + 2k_{3}(n) + k_{4}(n) )$ (1-5)

連立常微分方程式を考える

$\{\ \\ \frac{dy_1}{dx} = f_{1}( x,y_{1},y_{2},\cdots ,y_{n}) \\ \frac{dy_2}{dx} = f_{2}( x,y_{1},y_{2},\cdots ,y_{n}) \\ \frac{dy_3}{dx} = f_{3}( x,y_{1},y_{2},\cdots ,y_{n}) \\ \vdots \\ \frac{dy_n}{dx} = f_{n}( x,y_{1},y_{2},\cdots ,y_{n})$

ここでは、2元連立常微分方程式を解いてみる。

$\{\ \\ \frac{dy_1}{dx} = f_{1}( x,y_{1},y_{2}) \\ \frac{dy_2}{dx} = f_{2}( x,y_{1},y_{2})$

ここで、ルンゲクッタの公式(1-5)を適用すると、

$\LARGE y_{(1)n+1} = y_{(1)n} + \frac{1}{6}( k_{(1)1}(n) + 2k_{(1)2}(n) + 2k_{(1)3}(n) + k_{(1)4}(n) )$ (2-3)

$k_{(1)1}(n)=h*f_{1}(x_{n} , y_{(1)n} ,y_{(2)n} )$ (2-3a)
$k_{(1)2}(n)=h*f_{1}(x_n+\frac{h}{2}, y_{(1)n}+\frac{1}{2}*k_{(1)1}(n), y_{(2)n}+\frac{1}{2}*k_{(2)1}(n) )$ (2-3b)
$k_{(1)3}(n)=h*f_{1}(x_n+\frac{h}{2}, y_{(1)n}+\frac{1}{2}*k_{(1)2}(n) , y_{(2)n}+\frac{1}{2}*k_{(2)2}(n))$ (2-3c)
$k_{(1)4}(n)=h*f_{1}(x_n+h, y_{(1)n} + k_{(1)3}(n) , y_{(2)n} + k_{(2)3}(n) )$ (2-3d)

$\LARGE y_{(2)n+1} = y_{(2)n} + \frac{1}{6}( k_{(2)1}(n) + 2k_{(2)2}(n) + 2k_{(2)3}(n) + k_{(2)4}(n) )$ (2-4)

$k_{(2)1}(n)=h*f_{2}(x_{n} , y_{(1)n} ,y_{(2)n} )$ (2-4a)
$k_{(2)2}(n)=h*f_{2}(x_n+\frac{h}{2}, y_{(1)n}+\frac{1}{2}*k_{(1)1}(n), y_{(2)n}+\frac{1}{2}*k_{(2)1}(n) )$ (2-4b)
$k_{(2)3}(n)=h*f_{2}(x_n+\frac{h}{2}, y_{(1)n}+\frac{1}{2}*k_{(1)2}(n) , y_{(2)n}+\frac{1}{2}*k_{(2)2}(n))$ (2-4c)
$k_{(2)4}(n)=h*f_{2}(x_n+h, y_{(1)n} + k_{(1)3}(n) , y_{(2)n} + k_{(2)3}(n) )$ (2-4d)

[連立常微分方程式で4次のルンゲクッタを使用する例題]

プログラムコード[連立常微分方程式で4次のルンゲクッタを使用する例題]

#include "stdafx.h"
#include<iostream>
#include<fstream>
#include<io.h>
#include<math.h>
#include "ComplexC++.h"

#define PI 3.1415926 

const static double Em=1.0;                              //電圧源   

const static double R1=1.0;                               //抵抗
const static double R2=1.0;                               //抵抗

const static double C1=1.0;                               //コイル
const static double C2=1.0;                               //コイル

const static double freq=0.1;                            //周波数0.1Hz
const static double omega=2*PI*freq;                     //角周波数

const static double dh=(1/freq)/30.;                     //1周期の30分割

static double K3[3]={0,0,0};	
static double L3[3]={0,0,0};
static double K4[4]={0,0,0,0};	
static double L4[4]={0,0,0,0};	

void rungeKutta3ziK(double t,double V1,double V2);
void rungeKutta3ziL(double t,double V1,double V2);

void rungeKutta4ziK(double t,double V1,double V2);
void rungeKutta4ziL(double t,double V1,double V2);

double v1_func(double t,double V1,double V2);
double v2_func(double t,double V1,double V2);

///////// 理論計算 ///////////////
Complex c1(double t);
Complex c2(double t);
const static double lamda1=(-3-sqrt(5.))/2;
const static double lamda2=(-3+sqrt(5.))/2;
 
int _tmain(int argc, _TCHAR* argv[])
{
  
	errno_t err; 
	FILE *fp_output;
	
	//変数
	double t=0;          //実時間
	double v0=0;          //初期電圧
	double V1=0.0;
	double V2=0.0;
	double V111=0.0;
	double V222=0.0;

	Complex V11;
	Complex V22;
	

	if( err = fopen_s(&fp_output,"renritu_bibun_kekka.txt","w")  != 0){ exit(2);}
	
	t=0;
	for(int n=1; n< 3*(int)((1/freq)/dh); n++){

		//3次ルンゲクッタ係数の計算
		rungeKutta3ziK(t,V1,V2);
		rungeKutta3ziL(t,V1,V2);

		//4次ルンゲクッタ係数の計算
		rungeKutta4ziK(t,V111,V222);
		rungeKutta4ziL(t,V111,V222);
	
		//3次ルンゲクッタの計算
		V1=V1+(1./6.)*(K3[0]+4*K3[1]+K3[2])*dh;
		V2=V2+(1./6.)*(L3[0]+4*L3[1]+L3[2])*dh;

		//4次ルンゲクッタの計算
		V111=V111+(1./6.)*(K4[0]+2*K4[1]+2*K4[2]+K4[3]);
		V222=V222+(1./6.)*(L4[0]+2*L4[1]+2*L4[2]+L4[3]);
		
		//n=1から始まるから t=n*dh(n=1,2,3...)
		t=(double)n*dh;

		

		//電源電圧
		v0=1.+sin(omega*t);

		V11=c1(t)*2.*exp(lamda1*t)+c2(t)*2.*exp(lamda2*t);
		V22=c1(t)*(1-sqrt(5.))*exp(lamda1*t)+c2(t)*(1+sqrt(5.))*exp(lamda2*t);

		//printf("t = %f e=%f v1 = %f v2 = %f V11 = %f V22= %f\n",t,v0,V1,V2,cabs(V11),Real(V22));
		//fprintf_s(fp_output,"%f %f %f %f %f %f %f\n",t,v0,V1,V2,cabs(V11),Real(V22));

		printf("t = %f e=%f v1 = %f v2 = %f V111 = %f V222= %f V1real = %f V2real = %f\n",t,v0,V1,V2,V111,V222,Real(V11),Real(V22));
		fprintf_s(fp_output,"%f %f %f %f %f %f %f %f %f\n",t,v0,V1,V2,V111,V222,Real(V11),Real(V22));
	}

	if(fp_output != NULL) fclose(fp_output);
	getchar();
	return 0;
}   


void rungeKutta3ziK(double t,double V1,double V2){
	
	//これは t=n*dh(n=1,2,3,...)の計算となる
	K3[0]=v1_func(t,V1,V2);
	K3[1]=v1_func(t+dh/2,V1+K3[0]/2.0,V2+L3[0]/2.0);
	K3[2]=v1_func(t+dh,V1-K3[0]+2.0*K3[1],V2-L3[0]+2.0*L3[1]);
} 

void rungeKutta3ziL(double t,double V1,double V2){
	
	//これは t=n*dh(n=1,2,3,...)の計算となる
	L3[0]=v2_func(t,V1,V2);
	L3[1]=v2_func(t+dh/2,V1+K3[0]/2.0,V2+L3[0]/2.0);
	L3[2]=v2_func(t+dh,V1-K3[0]+2.0*K3[1],V2-L3[0]+2.0*L3[1]);
}

 
void rungeKutta4ziK(double t,double V1,double V2){
	//これは t=n*dh(n=1,2,3,...)の計算となる
 	K4[0]=dh*v1_func(t,V1,V2);
 	K4[1]=dh*v1_func(t+dh/2,V1+K4[0]/2.0,V2+L4[0]/2.0);
	K4[2]=dh*v1_func(t+dh/2,V1+K4[1]/2.0,V2+L4[1]/2.0);

	K4[3]=dh*v1_func(t+dh,V1+K4[2],V2+L4[2]);
}

void rungeKutta4ziL(double t,double V1,double V2){
	
 	//これは t=n*dh(n=1,2,3,...)の計算となる
	L4[0]=dh*v2_func(t,V1,V2);
	L4[1]=dh*v2_func(t+dh/2,V1+K4[0]/2.0,V2+L4[0]/2.0);
	L4[2]=dh*v2_func(t+dh/2,V1+K4[1]/2.0,V2+L4[1]/2.0);

	L4[3]=dh*v2_func(t+dh,V1+K4[2],V2+L4[2]);
}

Matlabでの描画

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 連立のルンゲクッタ法
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

load 'renritu_bibun_kekka.txt'

t=renritu_bibun_kekka(:,1);
e=renritu_bibun_kekka(:,2);
V1=renritu_bibun_kekka(:,3);
V2=renritu_bibun_kekka(:,4);
V111=renritu_bibun_kekka(:,5);
V222=renritu_bibun_kekka(:,6);
V11=renritu_bibun_kekka(:,7);
V22=renritu_bibun_kekka(:,8);

plot(t,e,'r-','linewidth',2);
hold on
plot(t,V1,'bo','linewidth',2);
hold on
plot(t,V111,'b+','linewidth',2);
hold on
plot(t,V11,'b-','linewidth',2);
hold on

plot(t,V2,'mo','linewidth',2);
hold on
plot(t,V222,'m+','linewidth',2);
hold on
plot(t,V22,'m-','linewidth',2);
hold on

grid on

xlabel('時間(sec)','Fontsize',18)
ylabel('電流(A)','Fontsize',18)

h=legend('電圧',...,
'V1(3次ルンゲクッタ)',...,
'V111(4次ルンゲクッタ)',...,
'V11(厳密解)',...,
'V2(3次ルンゲクッタ)',...,
'V222(4次ルンゲクッタ)',...,
'V22(厳密解)',...,
1);
set(h,'FontSize',15)

axis([0 Inf -Inf Inf]);
h=gca
get(h)
set(h,'LineWidth',2,...,
'FontSize',18)

print -djpeg renritu_method_byouga.jpg

3次のルンゲクッタ法

$\{\ \\ \frac{dy_1}{dx} = f_{1}( x,y_{1},y_{2}) \\ \frac{dy_2}{dx} = f_{2}( x,y_{1},y_{2})$

$\LARGE y_{(1)n+1} = y_{(1)n} + \frac{1}{6}( k_{(1)1}(n) + 4k_{(1)2}(n) + k_{(1)3}(n) )$ (3-1)

$k_{(1)1}(n)=h*f_{1}(x_{n} , y_{(1)n} ,y_{(2)n} )$ (3-1a)
$k_{(1)2}(n)=h*f_{1}(x_n+\frac{h}{2}, y_{(1)n}+\frac{1}{2}*k_{(1)1}(n), y_{(2)n}+\frac{1}{2}*k_{(2)1}(n) )$ (3-1b)
$k_{(1)3}(n)=h*f_{1}(x_n+ h , y_{(1)n} - k_{(1)1}(n) + 2k_{(1)2}(n) , y_{(2)n} - k_{(2)1}(n) + 2k_{(2)2}(n))$ (3-1c)

$\LARGE y_{(2)n+1} = y_{(2)n} + \frac{1}{6}( k_{(2)1}(n) + 4k_{(2)2}(n) + k_{(2)3}(n) )$ (3-2)

$k_{(2)1}(n)=h*f_{2}(x_{n} , y_{(1)n} ,y_{(2)n} )$ (3-2a)
$k_{(2)2}(n)=h*f_{2}(x_n+\frac{h}{2}, y_{(1)n}+\frac{1}{2}*k_{(1)1}(n), y_{(2)n}+\frac{1}{2}*k_{(2)1}(n) )$ (3-2b)
$k_{(2)3}(n)=h*f_{2}(x_n+ h , y_{(1)n} - k_{(1)1}(n) + 2k_{(1)2}(n) , y_{(2)n} - k_{(2)1}(n) + 2k_{(2)2}(n))$ (3-2c)