数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。noteにも書いています。https://note.mu/kumag_min

【連立一次方程式】双共役勾配法(BiCG法) C++コード付き

今回は正定値対称行列にしか適用できない共役勾配法(Conjugate Gradient法)をパワーアップした手法である双共役勾配法(Bi-Conjugate Gradient法)のC++コードを公開します。CG法についてはこちらをご覧下さい。

しくみはおいおい説明するとして、とりあえずコードをおいておきます。というより現在進行形で勉強中です!お待ちください。


今回解く連立一次方程式は

 Ax=b

とあらわされます。ここで、  A 3 \times 3 の行列、 x b 3 次元のベクトルです。 A b が既知で、 x が未知です。具体的には

 \left( \begin{array}{ccc} 3 & 1 & 1 \\ 1 & 4 & 1 \\ 1 & 1 & 5 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 5 \\ 6 \\ 7 \end{array} \right)

であり、これを解くと  x_1=1, x_2=1, x_3=1 となります。コードを書く際に、

数値計算のわざ

数値計算のわざ

の「第4章 大型線型方程式の反復解法」を参考しています。ではコードを以下に示します。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <stdio.h>
#include <time.h>
#include <stdlib.h> 

using namespace std;

int i,j,k,y,x;

//matix vector multiplication Ax=b
inline void mvm(double A[],double x[],double b[], int size)
{
	for(i=0;i<size;i++)
	{
		b[i]=0.;
		for(j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

//matix vector multiplication Atx=b (transpose)
inline void tmvm(double At[],double x[],double b[], int size)
{
	for(j=0;j<size;j++)
	{
		b[j]=0.;
		for(i=0;i<size;i++)
		{
			b[j]+=At[j*size+i]*x[i];
		}
	}
}

//inner product//
inline double ip(double a[], double b[], int size)
{
	double value=0.;
	for(i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

main()
{
	int size=3;
	double alpha,beta;
	double r0,rk,rk1;
	double eps=pow(2.0,-50);
	
	double *r=new double[size];
	double *rs=new double[size];
	double *p=new double[size];
	double *b=new double[size];
	double *ps=new double[size];
	double *Atps=new double[size];
	double *x=new double[size];
	double *Ax=new double[size];
	double *Ap=new double[size];
	double *A=new double[size*size];

	A[0*size+0]=3.;A[0*size+1]=1.;A[0*size+2]=1.;
	A[1*size+0]=1.;A[1*size+1]=4.;A[1*size+2]=1.;
	A[2*size+0]=1.;A[2*size+1]=1.;A[2*size+2]=5.;
	
	b[0]=5.;b[1]=6.;b[2]=7.;

	//IC//
	for(i=0;i<size;i++) x[i]=0.;
	
	mvm(A,x,Ax,size);
	for(i=0;i<size;i++) r[i]=b[i]-Ax[i];
	for(i=0;i<size;i++) p[i]=r[i];
	for(i=0;i<size;i++) ps[i]=r[i];
	for(i=0;i<size;i++) rs[i]=r[i];
	
	r0=sqrt(ip(r,r,size));
	
	//main loop//
	for(k=0;k<10000;k++)
	{
		mvm(A,p,Ap,size);
		alpha=ip(rs,r,size)/ip(ps,Ap,size);
		rk=ip(rs,r,size);
		for(i=0;i<size;i++) x[i]=x[i]+alpha*p[i];
		for(i=0;i<size;i++) r[i]=r[i]-alpha*Ap[i];
		
		rk1=sqrt(ip(r,r,size));
		if(rk1/r0<=eps) break;
		
		tmvm(A,ps,Atps,size);
		for(i=0;i<size;i++) rs[i]=rs[i]-alpha*Atps[i];
		beta=ip(rs,r,size)/rk;
		for(i=0;i<size;i++) p[i]=r[i]+beta*p[i];
		for(i=0;i<size;i++) ps[i]=rs[i]+beta*ps[i];
	}
	
	for(i=0;i<size;i++) cout<<x[i]<<endl;
	
	delete [] r,rs,p,b,ps,Atps,x,Ax,Ap,A;

    return 0;
}


Bi-CG法を安定化した、安定化双共役勾配法(Bi-conjugate gradient stabilized法)はこちらをご覧ください。