Ｒによる種々の方程式の解法
Solving nonlinear simultaneous equations by R langauge
作花 一志（京都情報大学院大学）
Kazushi Sakka（The Kyoto College of Graduate Studies for Informatics）
Abstract

R 言語は統計処理，データ解析のみならず数値計算に非常に便利である。この小文ではグラフの描き方から始め，線形・非線形方程式の解き方，陰関数で与えられた連立方程式，実解をもたない方程式さらに 1 階微分方程式までのグラフィカルに解く例を示す。
This paper shows the methods and examples of non-linear equations,spontaneous equations,
1st order differential equations.

■ まずグラフ描画

方程式を解くという学習は中学校以来学んでいるが，公式によって解ける方程式は 1 次方程式と 2次方程式くらいなもので，きわめて例外である。
f(x)=0 を解くとは x=f-1(0) を求めること，すなわち
y=f(x) と x 軸との交点を求めることである。
そこでまずはグラフを描いてみよう。y=sin(x) を－ 5 ≦ x ≦ 5 で描く最も簡単な方法は google でｙ＝ sin(x)を検索することであるが R プログラミ

ングのコードはcurve(sin(x),-5,5)だけでよく，一般に関数 f(x) を定義してグラフを描くには

図1

sin(x) 実線と cos(x) 破線

f= function(x) sin(x)curve(f(x),-5,5)

その他の色は“pink”，“skyblue”のようにすればよとすればよい。次の第 1 行は関数定義である。

い。線種を lty= で表す，1（実線）2（破線）3（点f= function(x) sin(x)g= function(x) cos(x)curve(f(x),-5,5)
curve(g(x),-5,5,col=2,add=T)
abline(v=0,col=8)
abline(h=0,col=8)= の代わりに <- を使うことも多い。
またここに別の関数 cos（x）のグラフを重ねて描くには add=T を付け加えるが，これは重ね描きを表す重要な命令で，付けないと前に描いたグラフは消えてしまう。col ＝の後の数はカラーコードで1（黒 : 通常省略），2（赤），3（緑），4（青），5（シアン），6（マゼンタ），7（黄），8（グレー）である。

48 NAIS Journal線）。また lwd= で線の幅を指定する。数値が大きいほど太い。なおパラメータ a,b を含めて
f=function (a,b,x) a*sin(x)+b*cos(x) と定義してcurve(f(1,0,x),-5,5)で sin(x) をcurve(f(0,1,x),-5,5，lty=2, add=T)
で cos(x) を描画してもよい。

■ 陰関数の場合■ 極座標の場合

関数を y=f(x) の形で与えられる場合を陽関数，グラフを表す式は絶対値 r と偏角 t を用いてf(x,y)=0 の形を陰関数という f(x,y)=0 より y=f(x) のr=f(t) と与えられ，閉曲線を描くとき便利である。

形に変換することはしばしば困難である。例えば

r=2*sin(2*t) の時

x +y =4 は容易にｙ＝±√ (4-x ) と変更できるが 2x=rcos(t)222

y=rsin(t) をプロットするには個の関数が必要であるし 3 次の項が含まれると大変だ。そこで陰関数の形のままグラフ表示してみよ

t=seq(0,2*pi,by=0.05)

う。 そ れ に は 等 高 線 を 描 く 関 数 contour を 使 う。

r=2*sin(2*t)

x -12xy -y +28=0（f1=0）と x -xy +y -12=0（f2=0x=r*cos(t);y=r*sin(t)楕円）を描くコードと結果グラフ図 2 を示す。

plot(x,y,xlim=c(-2,2),ylim=c(-2,2),type="l",xlab=

3322

ここで f2 として単にｙとすると f2=0 のグラフはy 軸と重なる"x",ylab="y")
これは正葉曲線といわれる。

######## 陰関数グラフ ############
x =y=seq(-10, 10, 0.2)
f1=function(x,y) x^3-12*x*y -y^3+28
f2=function(x,y) x^2-x*y +y^2-12
z1 = outer(x, y, f1)
z2 = outer(x, y, f2)
contour(x, y, z1, levels=0,col=2)
contour(x, y, z2, levels=0,col=4,add=T)
abline(h=0,col=8); abline(v=0,col=8)
grid()text(-5,9,"x^3-12*x*y -y^3+28=0",col=2)
text(5,9,"x^2-x*y +y^2-12=0",col=4)

図3正葉曲線■ 方程式解法
前 述 の と お り 方 程 式 f(x)=0 を 解 く と は 逆 変 換
x=f-1(0) を求めること，すなわち y=f(x) と x 軸との交点を求めることである
では実際に方程式を解いてみよう。
線形でも非線形でも同じことだから，いきなり非線形方程式 x2 － 2x=0 を解いてみよう。2x をｙと置き換えて・・・などとしてもうまくいかない。置換が悪いのではなく，そもそもいわゆる数学的解法
図2

陰関数表示が無理なのだ。受験技術に長けた大学生よりむしろ直観力のある小中学生のほうが早く解くことができるだろう。x=2 としてみると 22 － 22=0 だから解でz=locator(2)points(z,col=4)
p1=z$x[1] ; p2=z$x[2]ans=uniroot(f,c(p1,p2))# グラフと x 軸との交点を挟む 2 点をクリック# クリックした 2 点が表示される# クリックした 2 点の x 座標
# その間の解NAIS Journal49

あることはすぐにわかり，また x=4 も明らかに解

これにより関数を変更すればどんな方程式でも微である。しかし解はこの 2 つだけだろうか。これ分の知識なし解ける。ただしもちろん実根だけ。1以上は直観でも解析でもわからない。

点クリックで解く方法は後述する。

グラフを描いてみるともう一つ負の解があるが，これは数値的にしか求まらない。その解法は昔から多数開発されているが，R では uniroot という便利な関数が装備されている。x 軸との交点の左右に適整方程式の場合は係数をベクトルで与えてpolyroot 関数を使えばよい
x2+2x-3=0 の解はpolyroot(c(-3,2,1))

当な数 x1 と x2 を取り
ans = uniroot(f, c(x1, x2)) とするだけでよい。
解は ans$root に収まっていて，数値 -0.7666825より求まる。ただし（ ）内は昇べき順の係数である。
解は 1+0i-3+0iと求まるがこれは x=1，-3 を表す。

はコンソールとグラフィック画面に表示される。な
お四捨五入して小数点第 3 位まで表すなら

また 2x3+x2-2x+3=0 の解はround(ans$root,3) とする。

polyroot(c(3,-2,1,2))

より0.5767283+0.7580072i，# 方程式 f(x)=0の f(x) を与える

f = function(x) x^2-2^xcurve(f(x),-6,6)# この範囲の f(x) プロット

abline(h = 0, col = 8)

# x 軸を描く

abline(v = 0, col = 8) # ｙ 軸を描く
x1=;x2=# x1,x2 の値を与える

ans = uniroot(f, c(x1, x2)) # x1,x2 の間の解
ans$root#

text (-2,5,ans$root)

解はグラフィック画面に解はコンソールに

title(main="Solution x^2=2^x")-1.6534566+0.0000000i ，0.57672830.7580072iと求まるが 2 番目の値は実数で他は複素数（共役）である。
陰関数の場合は陽関数に変換せずに連立方程式の形で解を求めることができる。コードがやや長くなるがマウスクリックするのは 1 回だけでよい
x2-3x-3-y=0
xsin(-x)+2-y=0の解を求めよう
まず nleqslv というライブラリをインストールしておく。プログラムの第 9 行までは 2 つの曲線のと意義と描画を表している。次に関数 g を第 12 行～

でも x1,x2 にどんな値を与えたらいいか？グラフ

第 18 行のように定義して方程式を解く準備をする。

を描いてみてから x1=-2，x2=0 と x 軸との交点の左右の点を決めてもいいが，もっと素晴らしい方法

x =seq(-6, 6, 0.2);y=seq(-6, 25, 0.2)がある。locator を使えばグラフと x 軸との交点をf1=function(x,y) x^2-3*x-3-y

挟む 2 点をクリックするとその 2 点の x 座標がf2=function(x,y) x*sin(-x)+2-yx1,x2 となるのだ。

z1 = outer(x, y, f1)
z2 = outer(x, y, f2)
contour(x, y, z1, levels=0,col=2)
contour(x, y, z2, levels=0,col=4,add=T)
abline(h=0,col=8); abline(v=0,col=8)
grid()#define eqs
library(nleqslv)
g= function(z) {x = z[1]y = z[2]f1 = x^2-3*x-3-y
f2 = x*sin(-x)+2-y

図450 NAIS Journal

2xx =2 の解一 方 y=0 と し て 2 個 の 式 を ま と め てc(f1, f2)

x /2+x2-x-1=0 から前節のように 2 点クリックで求3}めたものが図 7 である。

#solvep=locator(1); points(p,col=4)x0=p$x ; y0=p$y
ans = nleqslv(c(x0,y0), g)
ans$x;points(ans$x[1],ans$x[2],pch=20,col=2)

当然解は同じである。
以上のことより陰関数で与えられた非線形連立方
程式でも実解を 1 点クリックでグラフィカルに求めることができる

arrows(x0,y0,ans$x[1],ans$x[2],lty=3)
text(x0,10,paste("x=",round(ans$x[1],3)))

図5

図7

1 点クリック

2 点クリック

曲線が描かれたグラフ上でクリックした点を(x0,y0) とすれば 2 曲線の交点が画面上に示される。
解は ans = nleqslv(c(x0,y0), g) より求まり
交点は x=ans$x[1]y=ans$x[2] である。

f1(x,y)=x3/2+x2-x(y+1)-1
f2(x,y)=yという連立方程式を上記の方法で解いた。f1=0は x=0 で不連続となり f2=0 はｙ軸と重なる

図81 点クリックによる連立方程式
x3-12xy -y3+28=0
x2-xy +y2-12=0の解は x=ans$x[1],y=ans$x[2]曲線は図 2 と同じ

図61 点クリック

NAIS Journal51

■ 導関数とニュートン法
関数を expression で定義し D を使うとその導関# ニュートン法の関数を定義
newton_method = function(fn1, fn2, x0) {eps=1e-7;max=100;x =x0数が得られる。複雑な公式を使わなくても微分でき

for (i in 1:max) {るので便利である。ただし関数形だけで関数値は求dx = fn1(x) / fn2(x)

まらないしグラフも描けない

x = x - dx
if (abs(dx) < eps) {y=expression (x^3+a*x+b*cos(x)+log(x)/a)(Dif1=D(y,"x"))#y(Dif2=D(Dif1,"x"))#y''(Dif3=D(Dif2,"x") )#y'''(Difa=D(y,"a"))#a で微分(Difb=D(y,"b"))#b で微分fn1(x) を微分してその導関数を求めるには
fn1=deriv(~ 関数形，“x”，func=T)
fn2(x) =function(x) attr(fn1(x),"gradient")とする。
fn2(x)=0 という方程式を解いて，その解に対する fn1(x) の値は極値である。
図 9 は fn1(x)-2sin(x)-x+2 の場合

return(x)}}stop("Not converge")}# 方程式とその導関数を定義
fn1=deriv(~x^2-2^x,"x",func=T)
fn2=function(x) attr(fn1(x),"gradient")
curve(fn1(x),-2,5)
abline(h=0,col=8) ;abline(v=0,col=8)# ニュートン法を適用
z=locator(1);points(z,col=4)x0=z$x(ans=newton_method(fn1, fn2, x0))
points(ans,0,col=2,pch=19)
title(paste("x=",round(ans,3)))

このプログラムは多少修正すれば実解を持たない方程式にも適応出来る。
例として x － log(x)-1=0 の場合，交点は存在しないがマウスクリックした点 z の x 座標，ｙ座標を実部虚部に持つ複素数をニュートン法の初期値とすればよい
x0=complex(re=z$x,imag=z$y)
これによってｚのｙ座標が正なら複素解
図9

導関数と極値

ニュートン法で方程式を解く場合には導関数が必要となるが関数形 fn1 を定義すればよい。次のプログラムでは⑥で関数を与え⑦で適当な点をクリックして初期値を与えて解いている。
なお最大反復回数を 100，収束判定定数を 10 － 7とした。

52 NAIS Journal-0.3951+1.7882iが求まり，負ならその共役の複 素 解 が 求 ま る。 同 様 に x2-x+2=0 の 解 は 0.5 ±1.3229i
であり，この値は当然 2 次方程式の解の公式から得られるものと一致する（図 10，図 11）。

なおライブラリ deSolve が必要である。
time1[492,] 49.1 3.141593[493,] 49.2 3.141593[494,] 49.3 3.141592[495,] 49.4 3.141593[496,] 49.5 3.141593[497,] 49.6 3.141593[498,] 49.7 3.141593[499,] 49.8 3.141593[500,] 49.9 3.141593[501,] 50.0 3.141593
図 10x-log（8x）+1=0の複素解
library(deSolve)#y’=x-y/2+2^(x)*sin(y)y0=1;fm="x-y/2+2^(x)*sin(y)"dfn=function(x,y,parms){list(x-y/2+2^(x)*sin(y))}
times = seq(from = 0, to = 50, by = 0.1)
out = ode(y = y0, times = times,
func = dfn,parms=null)
tail (out, n = 10)
plot(out[,1],out[,2],col=2,type="l",xlab="t",ylab="y")
図 11x2-x+2=0の複素解

title(paste("y'= ",fm,"y(0)=",y0))■ 一階微分方程式
微分方程式を解くとは x，y，y’，y’’ などを含む方程式から y を x の関数として表すことで，微分積分学修得の後で学ぶのが通例である。その起源はニュートンの運動方程式に始まり，これまでさまざまな解法が研究されているが，解析的方法は一般に非常に難解・技巧的であり，数値的にしか解けない場合も多い。
ところが R には ode と言う便利な関数が装備されていて，y’ と初期値を与えれば，非線形でも連立でも高階微分の場合でも容易に解

図 12

微分方程式の解初期値は（0,1）

くことができる。
例として非線形方程式 y’=x-y/2+2x*sin(y) を解いこの結果は初期値 y0=1 という特別な場合果なんてみよう。まず dfn で微分方程式 y’ を定義する。

だろうか？異なった初期値からでも同じ値が求まる

times では 0 から 50 まで 0.1 刻みで y，y’ の値をのか？異なった結果になるのか？それを調べるため

初 期 値 は x=0，y=1 と す る。 計 算 し out と い う

上記プログラムを一部変更して初期値は画面上を増

matrix に納める。out の第 1 列，第 2 列は t と y

すクリックして得ることを試みた。その結果初期値で最後の 10 行（tail）だけをコンソールに出力した。

y0>0 の場合はπに y0<0 の場合は - πに収束する

また plot によりグラフを描いた（図 12）。y の値はπに収束することがわかる。

ことが解った。図 13
このように微分方程式を短いコードで簡単に解くことができる。

NAIS Journal53

library(deSolve)plot(0,0,ylim=c(-3,3),type="n");abline(h=0,col=8)
z=locator(1)x0=z$x;y0=z$y

参考文献
https://www.f.waseda.jp/sakas/R/Rgraphics17.html
作花 NaisJ. Vol.15 p.74 2021
作花 NaisJ. Vol.17 p.41 2023

fm="x-y/2+2^(x)*sin(y)"dfn=function(x,y,parms){list(x-y/2+2^(x)*sin(y))}
times = seq(from = x0, to = 50, by = 0.1)
out = ode(y = y0, times = times, func = dfn,parms=null)
tail (out, n = 10)
plot(out[,1],out[,2],col=2,type="l",xlab="t",ylab="y")points(z,col=4,pch=20)
title(paste("y'= ",fm," ",out[nrow(out),2]))

図 13

54 NAIS Journal

微分方程式の解

初期値はマウスクリックで指定◆著者紹介

作花一志

Kazushi Sakka

京都情報大学院大学教授
京都大学理学研究科修了 理学博士
元京都コンピュータ学院鴨川校校長
元国際日本文化研究センター研究員
元日本天文教育普及研究会編集委員長
元京都大学総合人間学部非常勤講師
日本応用情報学会理事
暦学会理事