連立方程式

　　連立方程式 あるく 2007/11/21 18:11:37  (修正2回)
　　├ガウスの消去法は，通常，連立一次方程式を... 白石　和夫 2007/11/21 21:04:24 
　　│└いえ、ガウス法ではなく、ガウスの消去法の... あるく 2007/11/22 23:17:59 
　　│　└非線形の連立方程式を解くためのガウス消去... 白石　和夫 2007/11/23 09:06:48 
　　│　　└!２元非線形連立方程式の解 山中和義 2007/11/23 13:15:07  (修正1回)
　　│　　　├!オセッカイですが・・ SECOND 2007/11/24 00:21:32 
　　│　　　│└質問なのですが、皆さんが使ってる"!"はどう... あるく 2007/11/24 19:05:44 
　　│　　　│　└ヘルプファイルに、詳しい説明があります。 SECOND 2007/11/25 00:05:30  (修正1回)
　　│　　　├2.312376477871337.31237647787135 あるく 2007/11/27 17:43:16  (修正1回)
　　│　　　│└そうです。 山中和義 2007/11/27 18:44:39 
　　│　　　│　└ありがとうございます。疑問が解消しました... あるく 2007/11/27 19:17:46 
　　│　　　├自動微分でヤコビ行列を求める場合 山中和義 2007/11/28 21:30:22 
　　│　　　└たびたびすみません。また質問なのですが、 あるく 2007/12/02 16:39:20 
　　│　　　　└たとえば, 白石　和夫 2007/12/02 16:52:02  (修正1回)
　　│　　　　　└なるほど。 あるく 2007/12/02 18:08:53  (修正1回)
　　│　　　　　　└わかりやすい記述に変更してください。 山中和義 2007/12/03 10:04:15 
　　└!f(x)=x+5 SECOND 2007/11/22 18:47:13  (修正1回)
　　　└ありがとうございます。 あるく 2007/11/22 23:25:17 
　　　　└!複素数平面（ガウス平面）で根の見当をつけ... 山中和義 2007/11/23 13:52:37  (修正1回)

　　連立方程式 あるく 2007/11/21 18:11:37  (修正2回) ツリーへ

連立方程式 返事を書くノートメニュー

あるく <bydnumwjgh> 2007/11/21 18:11:37 ** この記事は2回修正されてます

f(x)=x+5
f(x)=√x

のような連立方程式を解くプログラムを作りたいのですが、ガウス法がどーのこーのでよくわかりません。

よかったら作り方を教えてください。

　　├ガウスの消去法は，通常，連立一次方程式を... 白石　和夫 2007/11/21 21:04:24  ツリーへ

Re: 連立方程式 返事を書くノートメニュー

白石　和夫 <ynwythjfwu> 2007/11/21 21:04:24

ガウスの消去法は，通常，連立一次方程式を解くための手法を指すものなので，それは使えません。「ガウス法」では，どのようにして解を求めるのですか？

　　│└いえ、ガウス法ではなく、ガウスの消去法の... あるく 2007/11/22 23:17:59  ツリーへ

Re: ガウスの消去法は，通常，連立一次方程式を... 返事を書くノートメニュー

あるく <bydnumwjgh> 2007/11/22 23:17:59

いえ、ガウス法ではなく、ガウスの消去法のことです。
すみません。

　　│　└非線形の連立方程式を解くためのガウス消去... 白石　和夫 2007/11/23 09:06:48  ツリーへ

Re: いえ、ガウス法ではなく、ガウスの消去法の... 返事を書くノートメニュー

白石　和夫 <ynwythjfwu> 2007/11/23 09:06:48

非線形の連立方程式を解くためのガウス消去法というのを私は知らないのですが，x+5=√xという方程式を解くのが目的であれば，
x=√x-5と変形して，適当な近似値から始めて x←SQR(x)-5 という代入計算を反復すれば収束して解が求まることもあります。
この方法は，必ず解が求まるとは限りませんし，解が求まる場合でもそのうちの一つしか求まりません。
10 OPTION ARITHMETIC COMPLEX
20 LET x=0
30 DO
40 LET x0=x
50 LET x=SQR(x)-5
60 ! PRINT x
70 LOOP UNTIL ABS((x-x0)/x)<1e-16
80 PRINT x
90 END

なお，複素数の範囲で考えるときには，√の意味を明確にしておく必要があります。十進BASICのSQR関数は，偏角が－πからπの間の平方根を表す約束になっていることに注意してください。

　　│　　└!２元非線形連立方程式の解 山中和義 2007/11/23 13:15:07  (修正1回) ツリーへ

Re: 非線形の連立方程式を解くためのガウス消去... 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/11/23 13:15:07 ** この記事は1回修正されてます

!２元非線形連立方程式の解

OPTION ARITHMETIC complex

DEF f(x,y)=x+5-y !f(x)=x+5
DEF g(x,y)=SQR(x)-y !g(x)=SQR(x)

LET d=0.1 !∂

LET xi=0 !初期値
LET yi=0

FOR i=1 TO 15 !収束させる

!連立１次方程式Jx=bを得る
DIM J(2,2)
LET fxy=f(xi,yi) !差分
LET J(1,1)=(f(xi+d,yi)-fxy)/d !J=(∂f/∂x　∂f/∂y)　ヤコビ行列
LET J(1,2)=(f(xi,yi+d)-fxy)/d !　(∂g/∂x　∂g/∂y)
LET gxy=g(xi,yi)
LET J(2,1)=(g(xi+d,yi)-gxy)/d
LET J(2,2)=(g(xi,yi+d)-gxy)/d

DIM b(2)
LET b(1)=-f(xi,yi) !b=(-f)
LET b(2)=-g(xi,yi) !　(-g)

!----------　※ガウスの消去法などで２元連立１次方程式を解く
DIM Ji(2,2) !逆行列を求める
MAT Ji=INV(J)

DIM x(2) !解⊿x、⊿y
MAT x=Ji*b
!----------

LET xi=xi+x(1) !Xi+1=Xi+⊿x
LET yi=yi+x(2) !Yi+1=Yi+⊿y

PRINT xi,yi !確認

NEXT i

END

たぶん、こんな感じでしょうか！？

　　│　　　├!オセッカイですが・・ SECOND 2007/11/24 00:21:32  ツリーへ

Re: !２元非線形連立方程式の解 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/11/24 00:21:32

!オセッカイですが・・
!このプログラムは、私を含め初等の学徒には、難しいと思いますので、翻訳して応援します。

!２元非線形連立方程式の解

!y=x+5
!y=SQR(x)

!0=x+5-y
!0=SQR(x)-y

!0=f(x,y)
!0=g(x,y)

!0=f(x,y), 0=g(x,y) の両方を満たすx,y が見つかれば解けたことになる。

!ここで、０になっていない X=f(x,y) と Y=g(x,y)で、抽象的な２次元ベクトル空間を考える。

!ベクトル(x,y) が、(⊿x,⊿y)だけ動いた時、ベクトル(X,Y)がどれ位動くかの微分係数のような
!ものがあれば、ニュートン法と同じ様に、

!出力ベクトル(X,Y)を、(0,0)にするための差分(0-X,0-Y) に対応する
!入力ベクトル(x,y)の、変化量(⊿x,⊿y)を逆算する事が、できる。

!ベクトルから、ベクトルへの微分係数は、行列の形式となります。（ヤコビ行列）

!　x,y に適当な初期値、０などをセットして、
!***** ここから *****

!入力列ベクトル(x,y) に対する出力列ベクトル( f(x,y),g(x,y) )への導関数(微分)としての行列
!（ヤコビ行列）で、一次の連立方程式を作り、入力差分(⊿x,⊿y)を逆算する。

!　|∂f(x,y)/∂x　∂f(x,y)/∂y| |⊿x|＝|0-f(x,y)|
!　|∂g(x,y)/∂x　∂g(x,y)/∂y| |⊿y|　|0-g(x,y)|

!両辺に、逆行列を乗ずると、左辺の行列は単位行列となって消える。(一次連立方程式の解法)

!　|⊿x|＝|∂f(x,y)/∂x　∂f(x,y)/∂y|(inv) |0-f(x,y)|
!　|⊿y|　|∂g(x,y)/∂x　∂g(x,y)/∂y|　　　|0-g(x,y)|

!　x=x+⊿x
!　y=y+⊿y

!***** ここまでを繰返す。*****

!そして、x,y が変化しなくなったら、そこが、出力ベクトル(X,Y)が、(0,0)に到着した所です。

!0=f(x,y)
!0=g(x,y)

!の条件を,満たしている入力ベクトル(x,y)です。

　　│　　　│└質問なのですが、皆さんが使ってる"!"はどう... あるく 2007/11/24 19:05:44  ツリーへ

Re: !オセッカイですが・・ 返事を書くノートメニュー

あるく <bydnumwjgh> 2007/11/24 19:05:44

質問なのですが、皆さんが使ってる"!"はどういう意味があるのですか？

　　│　　　│　└ヘルプファイルに、詳しい説明があります。 SECOND 2007/11/25 00:05:30  (修正1回) ツリーへ

Re: 質問なのですが、皆さんが使ってる"!"はどう... 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/11/25 00:05:30 ** この記事は1回修正されてます

ヘルプファイルに、詳しい説明があります。

目次→基本→行、下のほうの行末注釈の項。
目次→基本→REM文。

行末注釈記号"!"は、１文字で、
REM の代わりにも、なるため、多用しています。

　　│　　　├2.312376477871337.31237647787135 あるく 2007/11/27 17:43:16  (修正1回) ツリーへ

Re: !２元非線形連立方程式の解 返事を書くノートメニュー

あるく <bydnumwjgh> 2007/11/27 17:43:16 ** この記事は1回修正されてます

2.31237647787133 7.31237647787135
-6.27211452270137 -1.27211452270137
(-4.5665215844442 2.16253411536162) ( .433478415555792 2.16253411536161)
(-4.49999334095333 2.17942076048122) ( .500006659046668 2.17942076048122)
(-4.49999997350064 2.17944949328781) ( .500000026499363 2.17944949328781)
(-4.50000000003586 2.17944947178697) ( .499999999964139 2.17944947178697)
(-4.5 2.17944947177029) ( .500000000000001 2.17944947177029)
(-4.5 2.17944947177034) ( .5 2.17944947177034)
(-4.5 2.17944947177034) ( .5 2.17944947177034)
(-4.5 2.17944947177034) ( .5 2.17944947177034)
(-4.5 2.17944947177034) ( .5 2.17944947177034)
(-4.5 2.17944947177034) ( .5 2.17944947177034)
(-4.5 2.17944947177034) ( .5 2.17944947177034)
(-4.5 2.17944947177034) ( .5 2.17944947177034)
(-4.5 2.17944947177034) ( .5 2.17944947177034)

実行してみたらこのような値が出たのですが、左(例えば(-4.5 2.17944947177034))がxi、右(例えば( .5 2.17944947177034) )がyi
ということでいいんでしょうか？

　　│　　　│└そうです。 山中和義 2007/11/27 18:44:39  ツリーへ

Re: 2.312376477871337.31237647787135 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/11/27 18:44:39

そうです。

(-4.5 2.17944947177034)は、複素数の表示形式ですので、
実部が-4.5、虚部が2.17944947177034となります。

これは、SECONDさんの代数的解法の
　PRINT "根x=";(-9+SQR(9^2-4*1*25))/(2*1)
です。

実際に、解になっているかは、
　PRINT f(xi,yi),g(xi,yi) !０か、０に近ければ

をEND文の前に追加して確認できます。

　　│　　　│　└ありがとうございます。疑問が解消しました... あるく 2007/11/27 19:17:46  ツリーへ

Re: そうです。 返事を書くノートメニュー

あるく <bydnumwjgh> 2007/11/27 19:17:46

ありがとうございます。疑問が解消しました。

　　│　　　├自動微分でヤコビ行列を求める場合 山中和義 2007/11/28 21:30:22  ツリーへ

Re: !２元非線形連立方程式の解 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/11/28 21:30:22

自動微分でヤコビ行列を求める場合

精度の違いで多少収束が早い。

!２元非線形連立方程式の解

OPTION ARITHMETIC complex

!自動微分による微分係数の計算（ボトムアップ方式）
SUB DiffSQR(g(), f()) !合成 sqr(g)
LET f(1)=SQR(g(1)) !h(g)
FOR i=2 TO UBOUND(f)
LET f(i)=g(i)/(2*SQR(g(1))) !∂h(g)/∂g * ∂g/∂xi
NEXT i
END SUB

DIM x(3),y(3)
DATA 1,1,0 ![1, ∂f/∂x, ∂f/∂y]
MAT READ x
DATA 1,0,1 ![1, ∂g/∂x, ∂g/∂y]
MAT READ y

DIM c5(3) !定数5
DATA 5,0,0
MAT READ c5

FOR k=1 TO 15 !漸化式を収束させる

DIM xy(2) !⊿x、⊿y

!連立１次方程式Jx=bを得る
DIM v1(3),v2(3),v3(3),v4(3)
MAT v1=x+c5 !f(x,y)=x+5-y
MAT v2=v1-y
CALL DiffSQR(x, v3) !g(x,y)=SQR(x)-y
MAT v4=v3-y

DIM J(2,2)
LET J(1,1)=v2(2) !J=(∂f/∂x　∂f/∂y)　ヤコビ(Jacobi)行列
LET J(1,2)=v2(3) !　(∂g/∂x　∂g/∂y)
LET J(2,1)=v4(2)
LET J(2,2)=v4(3)

DIM b(2)
LET b(1)=-v2(1) !b=(-f)
LET b(2)=-v4(1) !　(-g)

!----------　※ガウスの消去法などで２元連立１次方程式を解く
DIM Ji(2,2) !逆行列を求める　Full BASIC版
MAT Ji=INV(J) !正則なら
MAT xy=Ji*b
!----------

LET x(1)=x(1)+xy(1) !Xi+1=Xi+⊿x
LET y(1)=y(1)+xy(2) !Yi+1=Yi+⊿y
PRINT x(1),y(1) !確認　∥Xi+1 - Xi∥≦ε(1+∥Xi+1∥)

NEXT k

DEF f(x,y)=x+5-y !f(x)=x+5
DEF g(x,y)=SQR(x)-y !g(x)=SQR(x)
PRINT f(x(1),y(1)),g(x(1),y(1)) !確認　∥f(xi)∥<δ

END

　　│　　　└たびたびすみません。また質問なのですが、 あるく 2007/12/02 16:39:20  ツリーへ

Re: !２元非線形連立方程式の解 返事を書くノートメニュー

あるく <bydnumwjgh> 2007/12/02 16:39:20

たびたびすみません。また質問なのですが、
DIM　の使い道がよく理解できません。このプログラムの中でしばしば見かけるのですが、どういった意図で使われているのですか?BASICヘルプには「添字の上限が～」みたいなことが書いてありましたが、添字って何？、という感じで理解できませんでした…

　　│　　　　└たとえば, 白石　和夫 2007/12/02 16:52:02  (修正1回) ツリーへ

Re: たびたびすみません。また質問なのですが、 返事を書くノートメニュー

白石　和夫 <ynwythjfwu> 2007/12/02 16:52:02 ** この記事は1回修正されてます

たとえば,
DIM a(2,3)
と書けば，aは2行3列の行列(matrix)になります。
aのi行j列成分をBASICではa(i,j)と書きます。
aの各構成要素a(1,1),a(1,2),a(1,3),a(2,1),a(2,2),a(2,3)は普通の変数と同様に使えますが，a(i,j)のi，jの部分に変数を含む数値式を書くことができるのが普通の変数との違いです。

　　│　　　　　└なるほど。 あるく 2007/12/02 18:08:53  (修正1回) ツリーへ

Re: たとえば, 返事を書くノートメニュー

あるく <bydnumwjgh> 2007/12/02 18:08:53 ** この記事は1回修正されてます

なるほど。
では
DIM Ji(2,2)
MAT Ji=INV(J)
この場合、Jiは2行2列の行列で、Jの逆行列である。という意味になるのでしょうか?
それと
DIM x(2)
この場合は2行1列になるのでしょうか?

　　│　　　　　　└わかりやすい記述に変更してください。 山中和義 2007/12/03 10:04:15  ツリーへ

Re: なるほど。 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/12/03 10:04:15

わかりやすい記述に変更してください。

>Jiは2行2列の行列で、Jの逆行列である。という意味になるのでしょうか?

そうです。

>DIM x(2)
>この場合は2行1列になるのでしょうか?

Ax=bから、これはベクトルですが、MAT文では、2行1列の行列として取り扱われます。
BASICでは、１次元の配列（添え字が１つ）です。

通常のプログラム言語では、このような「数の組」を配列を使って表現します。

下記のように、DIM x(2,1)と宣言したほうがいいかと思います。DIM b(2,1)も同様。

!２元非線形連立方程式の解

OPTION ARITHMETIC complex

DEF f(x,y)=x+5-y !f(x)=x+5
DEF g(x,y)=SQR(x)-y !g(x)=SQR(x)

LET d=0.1 !∂

LET xi=0 !初期値
LET yi=0

FOR i=1 TO 15 !収束させる

!連立１次方程式Jx=bを得る
DIM J(2,2)
LET fxy=f(xi,yi) !差分
LET J(1,1)=(f(xi+d,yi)-fxy)/d !J=(∂f/∂x　∂f/∂y)　ヤコビ行列
LET J(1,2)=(f(xi,yi+d)-fxy)/d !　(∂g/∂x　∂g/∂y)
LET gxy=g(xi,yi)
LET J(2,1)=(g(xi+d,yi)-gxy)/d
LET J(2,2)=(g(xi,yi+d)-gxy)/d

DIM b(2,1)
LET b(1,1)=-f(xi,yi) !b=(-f)
LET b(2,1)=-g(xi,yi) !　(-g)

!----------　※ガウスの消去法などで２元連立１次方程式を解く
DIM Ji(2,2) !逆行列を求める
MAT Ji=INV(J)

DIM x(2,1) !解⊿x、⊿y
MAT x=Ji*b
!----------

LET xi=xi+x(1,1) !Xi+1=Xi+⊿x
LET yi=yi+x(2,1) !Yi+1=Yi+⊿y

PRINT xi,yi !確認

NEXT i

END

　　└!f(x)=x+5 SECOND 2007/11/22 18:47:13  (修正1回) ツリーへ

Re: 連立方程式 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/11/22 18:47:13 ** この記事は1回修正されてます

!　f(x)=x+5
!　f(x)=√x

!２つのf(x)が、解を持つというのは、同じｘでの２つのf(x)が同値になる所、
!という事ですから、解法名は気にせず、グラフを書いてみましょう。

OPTION ARITHMETIC complex
SET WINDOW -10,10,-10,10
SET COLOR MIX(15) 0,0,1 !グリッドの色
DRAW grid

!　f(x)=x+5 のグラフ

FOR x=-10 TO +10 STEP 0.1
PLOT LINES:x, x+5; !実数部だけで、虚数部は区間内０、(描画は、省略されています。)
NEXT X
PLOT LINES

!　f(x)=√x のグラフ

SET LINE width 2 !太線
LET x=-10
LET xE=+10
DO
LET xB=x
LET reB=re(SQR(x))
LET imB=im(SQR(x))
LET x=x+0.1
SET LINE COLOR 1 !黒
PLOT LINES:xB,reB; x,re(SQR(x)) !実数部
SET LINE COLOR 4 !赤
PLOT LINES:xB,imB; x,im(SQR(x)) !虚数部
LOOP UNTIL xE<=x

!黒と赤の２本の線は、２つ合せて１つの値の線を意味します。

!２つのf(x)が、同値になるという事は、黒(実数部)と赤(虚数部)、
!両方が同時に一致する事です。
!ｘ＝－５で、黒(実数部)同士が、重なっているように見えますが、
!赤(虚数部)が、f(x)=x+5 の方は書かれていませんが、０で、同値でありません。

!x が実数であると、解が無いのが、わかります。
!代数的に、解いてみると、

!　x+5=√x
!　(x+5)^2=x　x^2+10x+25=x　x^2+9x+25=0
!　根 x＝ (-9±√(9^2-4*1*25) )/(2*1)　・・・(-9+i√19)/2 , (-9-i√19)/2

PRINT "根x =(実数虚数)=";(-9+SQR(9^2-4*1*25) )/(2*1)
PRINT "　　=(実数虚数)=";(-9-SQR(9^2-4*1*25) )/(2*1)

!やはり、複素数ｘが、解になっています。
!プログラムとしては、２次方程式の根の式を、置くぐらいしかないでしょう。

!これを、グラフ化するには、実数と虚数の２次元平面がｘ入力で、出力f(x)も、
!複素数ですので特別な工夫をしないと書けません。ご自身で工夫してみて下さい。

END

　　　└ありがとうございます。 あるく 2007/11/22 23:25:17  ツリーへ

Re: !f(x)=x+5 返事を書くノートメニュー

あるく <bydnumwjgh> 2007/11/22 23:25:17

ありがとうございます。
まだまだ初心者なので、完全にＳＥＣＯＮＤさんの書かれたことを理解できたわけではないのですが、がんばってみたいと思います。

　　　　└!複素数平面（ガウス平面）で根の見当をつけ... 山中和義 2007/11/23 13:52:37  (修正1回) ツリーへ

Re: ありがとうございます。 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/11/23 13:52:37 ** この記事は1回修正されてます

!複素数平面（ガウス平面）で根の見当をつける

OPTION ARITHMETIC complex

SET POINT STYLE 1

SET WINDOW -7,3,-5,5 !表示範囲
DRAW grid !座標

DEF f(x)=x+5-SQR(x) !x+5=SQR(x)
!DEF f(x)=(x+5)^2-x !x+5=SQR(x)の変形
!DEF f(x)=x^2+9*x+25 !x+5=SQR(x)の変形

FOR x=-7 TO 3 STEP 0.02 !複素数平面を走査する
FOR y=-5 TO 5 STEP 0.02
LET z=complex(x,y)
LET w=f(z)
IF ABS(1/w)>3 THEN PLOT POINTS: x,y !近似解
NEXT y
NEXT x

END

インデックスへ

EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。

連立方程式	返事を書くノートメニュー
あるく <bydnumwjgh> 2007/11/21 18:11:37 ** この記事は2回修正されてます
f(x)=x+5 f(x)=√x のような連立方程式を解くプログラムを作りたいのですが、ガウス法がどーのこーのでよくわかりません。よかったら作り方を教えてください。

Re: 連立方程式	返事を書くノートメニュー
白石　和夫 <ynwythjfwu> 2007/11/21 21:04:24
ガウスの消去法は，通常，連立一次方程式を解くための手法を指すものなので，それは使えません。「ガウス法」では，どのようにして解を求めるのですか？

Re: ガウスの消去法は，通常，連立一次方程式を...	返事を書くノートメニュー
あるく <bydnumwjgh> 2007/11/22 23:17:59
いえ、ガウス法ではなく、ガウスの消去法のことです。すみません。

Re: いえ、ガウス法ではなく、ガウスの消去法の...	返事を書くノートメニュー
白石　和夫 <ynwythjfwu> 2007/11/23 09:06:48
非線形の連立方程式を解くためのガウス消去法というのを私は知らないのですが，x+5=√xという方程式を解くのが目的であれば， x=√x-5と変形して，適当な近似値から始めて x←SQR(x)-5 という代入計算を反復すれば収束して解が求まることもあります。この方法は，必ず解が求まるとは限りませんし，解が求まる場合でもそのうちの一つしか求まりません。 10 OPTION ARITHMETIC COMPLEX 20 LET x=0 30 DO 40 LET x0=x 50 LET x=SQR(x)-5 60 ! PRINT x 70 LOOP UNTIL ABS((x-x0)/x)<1e-16 80 PRINT x 90 END なお，複素数の範囲で考えるときには，√の意味を明確にしておく必要があります。十進BASICのSQR関数は，偏角が－πからπの間の平方根を表す約束になっていることに注意してください。

Re: 非線形の連立方程式を解くためのガウス消去...	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/11/23 13:15:07 ** この記事は1回修正されてます
!２元非線形連立方程式の解 OPTION ARITHMETIC complex DEF f(x,y)=x+5-y !f(x)=x+5 DEF g(x,y)=SQR(x)-y !g(x)=SQR(x) LET d=0.1 !∂ LET xi=0 !初期値 LET yi=0 FOR i=1 TO 15 !収束させる !連立１次方程式Jx=bを得る DIM J(2,2) LET fxy=f(xi,yi) !差分 LET J(1,1)=(f(xi+d,yi)-fxy)/d !J=(∂f/∂x　∂f/∂y)　ヤコビ行列 LET J(1,2)=(f(xi,yi+d)-fxy)/d !　(∂g/∂x　∂g/∂y) LET gxy=g(xi,yi) LET J(2,1)=(g(xi+d,yi)-gxy)/d LET J(2,2)=(g(xi,yi+d)-gxy)/d DIM b(2) LET b(1)=-f(xi,yi) !b=(-f) LET b(2)=-g(xi,yi) !　(-g) !----------　※ガウスの消去法などで２元連立１次方程式を解く DIM Ji(2,2) !逆行列を求める MAT Ji=INV(J) DIM x(2) !解⊿x、⊿y MAT x=Ji*b !---------- LET xi=xi+x(1) !Xi+1=Xi+⊿x LET yi=yi+x(2) !Yi+1=Yi+⊿y PRINT xi,yi !確認 NEXT i END たぶん、こんな感じでしょうか！？

Re: !２元非線形連立方程式の解	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/11/24 00:21:32
!オセッカイですが・・ !このプログラムは、私を含め初等の学徒には、難しいと思いますので、翻訳して応援します。 !２元非線形連立方程式の解 !y=x+5 !y=SQR(x) !0=x+5-y !0=SQR(x)-y !0=f(x,y) !0=g(x,y) !0=f(x,y), 0=g(x,y) の両方を満たすx,y が見つかれば解けたことになる。 !ここで、０になっていない X=f(x,y) と Y=g(x,y)で、抽象的な２次元ベクトル空間を考える。 !ベクトル(x,y) が、(⊿x,⊿y)だけ動いた時、ベクトル(X,Y)がどれ位動くかの微分係数のような !ものがあれば、ニュートン法と同じ様に、 !出力ベクトル(X,Y)を、(0,0)にするための差分(0-X,0-Y) に対応する !入力ベクトル(x,y)の、変化量(⊿x,⊿y)を逆算する事が、できる。 !ベクトルから、ベクトルへの微分係数は、行列の形式となります。（ヤコビ行列） !　x,y に適当な初期値、０などをセットして、 !*** ここから * !入力列ベクトル(x,y) に対する出力列ベクトル( f(x,y),g(x,y) )への導関数(微分)としての行列 !（ヤコビ行列）で、一次の連立方程式を作り、入力差分(⊿x,⊿y)を逆算する。 !　\|∂f(x,y)/∂x　∂f(x,y)/∂y\| \|⊿x\|＝\|0-f(x,y)\| !　\|∂g(x,y)/∂x　∂g(x,y)/∂y\| \|⊿y\|　\|0-g(x,y)\| !両辺に、逆行列を乗ずると、左辺の行列は単位行列となって消える。(一次連立方程式の解法) !　\|⊿x\|＝\|∂f(x,y)/∂x　∂f(x,y)/∂y\|(inv) \|0-f(x,y)\| !　\|⊿y\|　\|∂g(x,y)/∂x　∂g(x,y)/∂y\|　　　\|0-g(x,y)\| !　x=x+⊿x !　y=y+⊿y !* ここまでを繰返す。*** !そして、x,y が変化しなくなったら、そこが、出力ベクトル(X,Y)が、(0,0)に到着した所です。 !0=f(x,y) !0=g(x,y) !の条件を,満たしている入力ベクトル(x,y)です。

Re: !オセッカイですが・・	返事を書くノートメニュー
あるく <bydnumwjgh> 2007/11/24 19:05:44
質問なのですが、皆さんが使ってる"!"はどういう意味があるのですか？

Re: 質問なのですが、皆さんが使ってる"!"はどう...	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/11/25 00:05:30 ** この記事は1回修正されてます
ヘルプファイルに、詳しい説明があります。目次→基本→行、下のほうの行末注釈の項。目次→基本→REM文。行末注釈記号"!"は、１文字で、 REM の代わりにも、なるため、多用しています。

Re: !２元非線形連立方程式の解	返事を書くノートメニュー
あるく <bydnumwjgh> 2007/11/27 17:43:16 ** この記事は1回修正されてます
2.31237647787133 7.31237647787135 -6.27211452270137 -1.27211452270137 (-4.5665215844442 2.16253411536162) ( .433478415555792 2.16253411536161) (-4.49999334095333 2.17942076048122) ( .500006659046668 2.17942076048122) (-4.49999997350064 2.17944949328781) ( .500000026499363 2.17944949328781) (-4.50000000003586 2.17944947178697) ( .499999999964139 2.17944947178697) (-4.5 2.17944947177029) ( .500000000000001 2.17944947177029) (-4.5 2.17944947177034) ( .5 2.17944947177034) (-4.5 2.17944947177034) ( .5 2.17944947177034) (-4.5 2.17944947177034) ( .5 2.17944947177034) (-4.5 2.17944947177034) ( .5 2.17944947177034) (-4.5 2.17944947177034) ( .5 2.17944947177034) (-4.5 2.17944947177034) ( .5 2.17944947177034) (-4.5 2.17944947177034) ( .5 2.17944947177034) (-4.5 2.17944947177034) ( .5 2.17944947177034) 実行してみたらこのような値が出たのですが、左(例えば(-4.5 2.17944947177034))がxi、右(例えば( .5 2.17944947177034) )がyi ということでいいんでしょうか？

Re: 2.312376477871337.31237647787135	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/11/27 18:44:39
そうです。 (-4.5 2.17944947177034)は、複素数の表示形式ですので、実部が-4.5、虚部が2.17944947177034となります。これは、SECONDさんの代数的解法の　PRINT "根x=";(-9+SQR(9^2-4125))/(2*1) です。実際に、解になっているかは、　PRINT f(xi,yi),g(xi,yi) !０か、０に近ければをEND文の前に追加して確認できます。

Re: そうです。	返事を書くノートメニュー
あるく <bydnumwjgh> 2007/11/27 19:17:46
ありがとうございます。疑問が解消しました。

Re: !２元非線形連立方程式の解	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/11/28 21:30:22
自動微分でヤコビ行列を求める場合精度の違いで多少収束が早い。 !２元非線形連立方程式の解 OPTION ARITHMETIC complex !自動微分による微分係数の計算（ボトムアップ方式） SUB DiffSQR(g(), f()) !合成 sqr(g) LET f(1)=SQR(g(1)) !h(g) FOR i=2 TO UBOUND(f) LET f(i)=g(i)/(2SQR(g(1))) !∂h(g)/∂g ∂g/∂xi NEXT i END SUB DIM x(3),y(3) DATA 1,1,0 ![1, ∂f/∂x, ∂f/∂y] MAT READ x DATA 1,0,1 ![1, ∂g/∂x, ∂g/∂y] MAT READ y DIM c5(3) !定数5 DATA 5,0,0 MAT READ c5 FOR k=1 TO 15 !漸化式を収束させる DIM xy(2) !⊿x、⊿y !連立１次方程式Jx=bを得る DIM v1(3),v2(3),v3(3),v4(3) MAT v1=x+c5 !f(x,y)=x+5-y MAT v2=v1-y CALL DiffSQR(x, v3) !g(x,y)=SQR(x)-y MAT v4=v3-y DIM J(2,2) LET J(1,1)=v2(2) !J=(∂f/∂x　∂f/∂y)　ヤコビ(Jacobi)行列 LET J(1,2)=v2(3) !　(∂g/∂x　∂g/∂y) LET J(2,1)=v4(2) LET J(2,2)=v4(3) DIM b(2) LET b(1)=-v2(1) !b=(-f) LET b(2)=-v4(1) !　(-g) !----------　※ガウスの消去法などで２元連立１次方程式を解く DIM Ji(2,2) !逆行列を求める　Full BASIC版 MAT Ji=INV(J) !正則なら MAT xy=Ji*b !---------- LET x(1)=x(1)+xy(1) !Xi+1=Xi+⊿x LET y(1)=y(1)+xy(2) !Yi+1=Yi+⊿y PRINT x(1),y(1) !確認　∥Xi+1 - Xi∥≦ε(1+∥Xi+1∥) NEXT k DEF f(x,y)=x+5-y !f(x)=x+5 DEF g(x,y)=SQR(x)-y !g(x)=SQR(x) PRINT f(x(1),y(1)),g(x(1),y(1)) !確認　∥f(xi)∥<δ END

Re: !２元非線形連立方程式の解	返事を書くノートメニュー
あるく <bydnumwjgh> 2007/12/02 16:39:20
たびたびすみません。また質問なのですが、 DIM　の使い道がよく理解できません。このプログラムの中でしばしば見かけるのですが、どういった意図で使われているのですか?BASICヘルプには「添字の上限が～」みたいなことが書いてありましたが、添字って何？、という感じで理解できませんでした…

Re: たびたびすみません。また質問なのですが、	返事を書くノートメニュー
白石　和夫 <ynwythjfwu> 2007/12/02 16:52:02 ** この記事は1回修正されてます
たとえば, DIM a(2,3) と書けば，aは2行3列の行列(matrix)になります。 aのi行j列成分をBASICではa(i,j)と書きます。 aの各構成要素a(1,1),a(1,2),a(1,3),a(2,1),a(2,2),a(2,3)は普通の変数と同様に使えますが，a(i,j)のi，jの部分に変数を含む数値式を書くことができるのが普通の変数との違いです。

Re: たとえば,	返事を書くノートメニュー
あるく <bydnumwjgh> 2007/12/02 18:08:53 ** この記事は1回修正されてます
なるほど。では DIM Ji(2,2) MAT Ji=INV(J) この場合、Jiは2行2列の行列で、Jの逆行列である。という意味になるのでしょうか? それと DIM x(2) この場合は2行1列になるのでしょうか?

Re: なるほど。	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/12/03 10:04:15
わかりやすい記述に変更してください。 >Jiは2行2列の行列で、Jの逆行列である。という意味になるのでしょうか? そうです。 >DIM x(2) >この場合は2行1列になるのでしょうか? Ax=bから、これはベクトルですが、MAT文では、2行1列の行列として取り扱われます。 BASICでは、１次元の配列（添え字が１つ）です。通常のプログラム言語では、このような「数の組」を配列を使って表現します。下記のように、DIM x(2,1)と宣言したほうがいいかと思います。DIM b(2,1)も同様。 !２元非線形連立方程式の解 OPTION ARITHMETIC complex DEF f(x,y)=x+5-y !f(x)=x+5 DEF g(x,y)=SQR(x)-y !g(x)=SQR(x) LET d=0.1 !∂ LET xi=0 !初期値 LET yi=0 FOR i=1 TO 15 !収束させる !連立１次方程式Jx=bを得る DIM J(2,2) LET fxy=f(xi,yi) !差分 LET J(1,1)=(f(xi+d,yi)-fxy)/d !J=(∂f/∂x　∂f/∂y)　ヤコビ行列 LET J(1,2)=(f(xi,yi+d)-fxy)/d !　(∂g/∂x　∂g/∂y) LET gxy=g(xi,yi) LET J(2,1)=(g(xi+d,yi)-gxy)/d LET J(2,2)=(g(xi,yi+d)-gxy)/d DIM b(2,1) LET b(1,1)=-f(xi,yi) !b=(-f) LET b(2,1)=-g(xi,yi) !　(-g) !----------　※ガウスの消去法などで２元連立１次方程式を解く DIM Ji(2,2) !逆行列を求める MAT Ji=INV(J) DIM x(2,1) !解⊿x、⊿y MAT x=Ji*b !---------- LET xi=xi+x(1,1) !Xi+1=Xi+⊿x LET yi=yi+x(2,1) !Yi+1=Yi+⊿y PRINT xi,yi !確認 NEXT i END

Re: 連立方程式	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/11/22 18:47:13 ** この記事は1回修正されてます
!　f(x)=x+5 !　f(x)=√x !２つのf(x)が、解を持つというのは、同じｘでの２つのf(x)が同値になる所、 !という事ですから、解法名は気にせず、グラフを書いてみましょう。 OPTION ARITHMETIC complex SET WINDOW -10,10,-10,10 SET COLOR MIX(15) 0,0,1 !グリッドの色 DRAW grid !　f(x)=x+5 のグラフ FOR x=-10 TO +10 STEP 0.1 PLOT LINES:x, x+5; !実数部だけで、虚数部は区間内０、(描画は、省略されています。) NEXT X PLOT LINES !　f(x)=√x のグラフ SET LINE width 2 !太線 LET x=-10 LET xE=+10 DO LET xB=x LET reB=re(SQR(x)) LET imB=im(SQR(x)) LET x=x+0.1 SET LINE COLOR 1 !黒 PLOT LINES:xB,reB; x,re(SQR(x)) !実数部 SET LINE COLOR 4 !赤 PLOT LINES:xB,imB; x,im(SQR(x)) !虚数部 LOOP UNTIL xE<=x !黒と赤の２本の線は、２つ合せて１つの値の線を意味します。 !２つのf(x)が、同値になるという事は、黒(実数部)と赤(虚数部)、 !両方が同時に一致する事です。 !ｘ＝－５で、黒(実数部)同士が、重なっているように見えますが、 !赤(虚数部)が、f(x)=x+5 の方は書かれていませんが、０で、同値でありません。 !x が実数であると、解が無いのが、わかります。 !代数的に、解いてみると、 !　x+5=√x !　(x+5)^2=x　x^2+10x+25=x　x^2+9x+25=0 !　根 x＝ (-9±√(9^2-4125) )/(21)　・・・(-9+i√19)/2 , (-9-i√19)/2 PRINT "根x =(実数虚数)=";(-9+SQR(9^2-4125) )/(21) PRINT "　　=(実数虚数)=";(-9-SQR(9^2-4125) )/(2*1) !やはり、複素数ｘが、解になっています。 !プログラムとしては、２次方程式の根の式を、置くぐらいしかないでしょう。 !これを、グラフ化するには、実数と虚数の２次元平面がｘ入力で、出力f(x)も、 !複素数ですので特別な工夫をしないと書けません。ご自身で工夫してみて下さい。 END

Re: !f(x)=x+5	返事を書くノートメニュー
あるく <bydnumwjgh> 2007/11/22 23:25:17
ありがとうございます。まだまだ初心者なので、完全にＳＥＣＯＮＤさんの書かれたことを理解できたわけではないのですが、がんばってみたいと思います。

Re: ありがとうございます。	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/11/23 13:52:37 ** この記事は1回修正されてます
!複素数平面（ガウス平面）で根の見当をつける OPTION ARITHMETIC complex SET POINT STYLE 1 SET WINDOW -7,3,-5,5 !表示範囲 DRAW grid !座標 DEF f(x)=x+5-SQR(x) !x+5=SQR(x) !DEF f(x)=(x+5)^2-x !x+5=SQR(x)の変形 !DEF f(x)=x^2+9*x+25 !x+5=SQR(x)の変形 FOR x=-7 TO 3 STEP 0.02 !複素数平面を走査する FOR y=-5 TO 5 STEP 0.02 LET z=complex(x,y) LET w=f(z) IF ABS(1/w)>3 THEN PLOT POINTS: x,y !近似解 NEXT y NEXT x END