HSPポータル
サイトマップ お問い合わせ


HSPTV!掲示板


未解決 解決 停止 削除要請

2016
0330
いろは行列による連立方程式2解決


いろは

リンク

2016/3/30(Wed) 12:27:51|NO.75119

今作っているソフトの一環で連立n元一次方程式を解くプログラムを作っています。
一応書いたには書いたのですが、どこかが間違っているようで上手く動きません。どこが間違っているか教えていただけないでしょうか。

#include "hspmath.as" title "連立n元一次方程式を解く" font "MS 明朝",15 pos 10,0:input n,25,25:objsize 50,20:pos 35,0:button "決定",*gyouretu stop *gyouretu if (clear):clrobj objID(0),objID(1):objsize 50,20:pos 35,0:button "決定",*gyouretu //オブジェクトの消去と描写 nn = n color 255,255,255:boxf:color:pos 10,25:mes ""+nn+"連立元一次方程式" objID(0) = stat objsize 25,25 repeat nn //係数行列 cnt1 = cnt cnt_1 =cnt*nn repeat nn cnt_2 = cnt a(cnt_1+cnt_2) = cnt_1+cnt_2 pos 50+cnt_2*25,50+cnt1*25 input a(cnt_1+cnt_2),25,25 loop loop repeat nn //定数項 b(cnt) = cnt pos 50+nn*25+25,50+cnt*25 input b(cnt),25,25 loop pos 50,50+nn*25:mes "係数行列":pos 50+nn*25+25,50+nn*25:mes "定数項" objsize 50,20:pos 20,50+nn*25+25:button "結果",*keka objID(1) = stat clear = 1 stop *keka color 255,255,255:boxf 0,nn*25+25,640,nn*25+25 //解の初期化 //det(A)を求める repeat nn cnt_3 = cnt repeat nn cnt_4 = cnt if (cnt_4+cnt_3 != nn)&(p1 != 1):dapd(cnt_4) = a(cnt_4*nn + cnt_4+cnt_3) if (cnt_4+cnt_3 = nn)|(p1 = 1):dapd(cnt_4) = a(cnt_4*nn + cnt_4+cnt_3-nn):p1=1 loop p1 = 0 repeat nn-1 if (cnt = 0):dap_d(cnt_3) = dapd(cnt) if (cnt = 0):dap_(cnt_3) = dap_d(cnt_3) * dapd(cnt+1) if (cnt != 0):dap_(cnt_3) = dap_(cnt_3) * dapd(cnt+1) loop loop repeat nn dap += dap_(cnt) loop repeat nn cnt_5 = cnt repeat nn cnt_6 = cnt if (cnt_6+cnt_5 != nn)&(m1 != 1):damd(cnt_6) = a(cnt_6*nn + nn-1-cnt_6-cnt_5) if (cnt_6+cnt_5 = nn)|(m1 = 1):damd(cnt_6) = a(cnt_6*nn + nn-1-cnt_6-cnt_5+nn):m1=1 loop m1 = 0 repeat nn-1 if (cnt = 0):dam_d(cnt_5) = damd(cnt) if (cnt = 0):dam_(cnt_5) = dam_d(cnt_5) * damd(cnt+1) if (cnt != 0):dam_(cnt_5) = dam_(cnt_5) * damd(cnt+1) loop loop repeat nn dam += dam_(cnt) loop deta = dap - dam //det(A)が求まる if (deta = 0)::dialog "det(A)=0のため、解がない、もしくは解が複数存在します",1,"解なし":stop //余因子行列を求める repeat nn*nn cnt_12 = cnt //行列の要素の定義 repeat nn-1 cnt11 = cnt*(nn-1) cnt_11 = cnt*nn if (cnt_12-cnt_11 <= nn-1)&(cnt_12-cnt_11 >= 0):r=1 repeat nn-1 if (fmod(cnt_12,nn) = fmod(cnt,nn))|(g = 1):a_(cnt11+cnt) = a(cnt_11+cnt+1):g=1 if (r = 1):a_(cnt11+cnt) = a(cnt_11+cnt+nn) if (g = 1)&(r = 1):a_(cnt11+cnt) = a(cnt_11+cnt+nn+1) if (g != 1)&(r != 1):a_(cnt11+cnt) = a(cnt_11+cnt) loop g = 0 loop r = 0 repeat nn-1 cnt_7 = cnt repeat nn-1 cnt_8 = cnt if (cnt_8+cnt_7 != nn-1)&(p2 != 1):ydapd(cnt_8) = a_(cnt_8*(nn-1) + cnt_8+cnt_7) if (cnt_8+cnt_7 = nn-1)|(p2 = 1):ydapd(cnt_8) = a_(cnt_8*(nn-1) + cnt_8+cnt_7-(nn-1)):p2=1 loop p2 = 0 repeat nn-1-1 if (cnt = 0):ydap_d(cnt_7) = ydapd(cnt) if (cnt = 0):ydap_(cnt_7) = ydap_d(cnt_7) * ydapd(cnt+1) if (cnt != 0):ydap_(cnt_7) = ydap_(cnt_7) * ydapd(cnt+1) loop loop repeat nn-1 ydap(cnt_12) += ydap_(cnt) loop repeat nn-1 cnt_9 = cnt repeat nn-1 cnt_10 = cnt if (cnt_10+cnt_9 != nn-1)&(m2 != 1):ydamd(cnt_10) = a_(cnt_10*(nn-1) + (nn-1)-1-cnt_10-cnt_9) if (cnt_10+cnt_9 = nn-1)|(m2 = 1):ydamd(cnt_10) = a_(cnt_10*(nn-1) + (nn-1)-1-cnt_10-cnt_9+(nn-1)):m2=1 loop m2 = 0 repeat nn-1-1 if (cnt = 0):ydam_d(cnt_9) = ydamd(cnt) if (cnt = 0):ydam_(cnt_9) = ydam_d(cnt_9) * ydamd(cnt+1) if (cnt != 0):ydam_(cnt_9) = ydam_(cnt_9) * ydamd(cnt+1) loop loop repeat nn-1 ydam(cnt_12) += ydam_(cnt) loop ydeta(cnt_12) = ydap(cnt_12)-ydam(cnt_12) //余因子行列が求まる(順番はまだ) loop repeat nn //余因子行列の順番も決定 cnt_13 = cnt*nn cnt13 = cnt repeat nn if (fmod(cnt,2) = 0)&(fmod(cnt13,2) = 0){ deta_(cnt_13+cnt) = ydeta(cnt*nn + cnt13) }else:deta_(cnt_13+cnt) = -ydeta(cnt*nn + cnt13) if (fmod(cnt,2) = 0)&(fmod(cnt13,2) != 0){ deta_(cnt_13+cnt) = -ydeta(cnt*nn + cnt13) }else:deta_(cnt_13+cnt) = ydeta(cnt*nn + cnt13) loop loop //逆行列算出 repeat nn*nn gyaku(cnt) = double(deta_(cnt))/double(deta) loop //解の算出 repeat nn cnt_14 = cnt*nn repeat nn x_(cnt+cnt_14) = double(gyaku(cnt+cnt_14)) * double(b(cnt)) loop loop repeat nn cnt_15 = cnt repeat nn x(cnt_15) += double(x_(cnt)) loop loop repeat nn //解の表示 color:mes "x"+cnt+"="+x(cnt)+" loop stop
こんな感じで逆行列から解を求める形にしたのですがどこを直せばいいのかわかりません。お願いします。



この記事に返信する


cats

リンク

2016/3/30(Wed) 19:33:03|NO.75127

変数名やアルゴリズムがよく分からないので
コードはちゃんと読めていませんが、いくつかアドバイスを。
まず、見たところ1次元配列を使っているようですが、
2次元配列を使うほうが絶対に楽です。
行列式を求めるところで変な結果が出たので、
行列式の計算から間違っていると思います。
行列式を求めるのに余因子行列を使って計算するのもいいですが、
スタックも再帰も無いHSPでは、ちょっとややこしいと思います。
一度行列を上三角行列に変換して対角成分を掛ければ行列式になります。

n = length(A) // 行列Aを上三角行列に変換 ratio = 0.0 i = 0 : j = 0 : k = 0 repeat n i = cnt repeat n j = cnt if (j > i) { ratio = A(j, i) / A(i, i) repeat n k = cnt A(j, k) -= ratio * A(i, k) loop } loop loop // 行列式を求める det = 1.0 repeat n det *= A(cnt, cnt) loop mes det stop
※上のサンプルを動かすにはdouble型のn次正方行列Aを用意してください。



いろは

リンク

2016/3/31(Thu) 23:13:28|NO.75141

無事に正常に動作をさせることに成功しました。ありがとうございました。



ONION software Copyright 1997-2023(c) All rights reserved.