上図は、曲線に対する、通過点指定の接線を求める場合の計算方法です。
接線を求める範囲の最初の点を初期位置とします、初期位置から、小さな間隔で、xの値を等間隔で大きくします、その時の通過指定点と、xの値による曲線上の位置の角度を求めます。
その角度の変化は、接戦の位置を境に、変化の方向が逆になります。
その逆になる点を検出し、逆になったらピッチを二つ戻し、ピッチを小さくして、再度計算をします。
ピッチが、指定した値より小さくなったらその位置を非線形方程式の曲線上の接線の位置とします。
非線形方程式の解放続き 指定点を通過する接線
非線形方程式による曲線への指定点通過の接線で指定点は曲線上に無い場合です。
上図は、曲線に対する、通過点指定の接線を求める場合の計算方法です。
接線を求める範囲の最初の点を初期位置とします、初期位置から、小さな間隔で、xの値を等間隔で大きくします、その時の通過指定点と、xの値による曲線上の位置の角度を求めます。
その角度の変化は、接戦の位置を境に、変化の方向が逆になります。
その逆になる点を検出し、逆になったらピッチを二つ戻し、ピッチを小さくして、再度計算をします。
ピッチが、指定した値より小さくなったらその位置を非線形方程式の曲線上の接線の位置とします。
プログラム
プログラムは、二次方程式と三次方程式の計算例となっています。
微分して、直線との連立方程式を解けば、接線の値を求めることが出来ますが、あくまでも計算例です。
プログラム中の方程式を入れ替えれば、任意の方程式とすることが出来ます。
初期位置aは、接線を検出する最初の位置です。
初期間隔pは、次の位置を計算間隔ですので次の位置はx=x+pとなります。
検索上限は、接線が無い場合、あるいは有っても非常に大きな値となる場合もあるので、その場合、計算を打ち切る最大のx値です。
初期間隔Pの値は、大きく過ぎると、角度を計算するのに、象限を利用するのですが、二象限に跨ってしまうことがあり、角度の方向を正しく計算できないことがあります。
小さくしすぎると、計算数が増えるので、適当に調整します。
グラフ最小値、最大値は、初期位置aあるいは、検束上限の値より大きな範囲に設定すると、その値でグラフが表示されます。
二次方程式を選んだ場合は、接線は通過指定点の位置によって、接線は二本になるか、無いとなります、三次方程式の場合は、一本か三本となります。
計算は、最初に求められた一本しか表示されません、次の接線を求める場合は、初期位置aの値を変更します。
又、二次方程式を選んだ場合だけ、連立方程式を解いた、xの値がt1,t2として表示されます、近似値との誤差を確認する為のものです。
近似値の精度は、doubleだとすると、その精度桁数の二分の一ぐらいです。
// 非線形方程式と通過点指定の接線計算 // 曲線上の座標をx方向へ移動し曲線上の点と指定通過点中心角度の回転方向から接線の傾斜の近似値を計算します。 // 計算例は、容易に微分出来る方程式です、あくまでも計算例です。 unit Main; interface uses Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics, Vcl.Controls, Vcl.Forms, Vcl.Dialogs, VclTee.TeeGDIPlus, VCLTee.Series, Vcl.StdCtrls, Vcl.Buttons, Vcl.ExtCtrls, VCLTee.TeEngine, VCLTee.TeeProcs, VCLTee.Chart; type TForm1 = class(TForm) Chart1: TChart; Series1: TLineSeries; LabeledEdit1: TLabeledEdit; LabeledEdit2: TLabeledEdit; LabeledEdit3: TLabeledEdit; LabeledEdit4: TLabeledEdit; LabeledEdit5: TLabeledEdit; LabeledEdit6: TLabeledEdit; LabeledEdit7: TLabeledEdit; BitBtn1: TBitBtn; Memo1: TMemo; Series3: TPointSeries; Series4: TPointSeries; BitBtn2: TBitBtn; Series2: TLineSeries; CheckBox1: TCheckBox; RadioButton1: TRadioButton; RadioButton2: TRadioButton; procedure BitBtn1Click(Sender: TObject); procedure BitBtn2Click(Sender: TObject); procedure CheckBox1Click(Sender: TObject); procedure FormCreate(Sender: TObject); procedure LabeledEditChange(Sender: TObject); procedure RadioButton1Click(Sender: TObject); procedure RadioButton2Click(Sender: TObject); private { Private 宣言 } procedure chart_graph(gl, gr: double); public { Public 宣言 } LabeledEdit : array[1..7] of TLabeledEdit; end; var Form1: TForm1; implementation {$R *.dfm} uses system.Math; // ****** 参考 計算による接線の値 ******************************************** // f(x)=(x-3)^2-2=x^2-6x+9-2=x^2-6x+7 // 微分 傾斜α // f'(x)=2x-6 // f(x)上の接点x座標をtとします。 // 指定点座標yの値は接点(x)t座標と指定点x座標の差分に傾斜αを掛けたものに // 接点座標yの値を加算したものになります。 // y=(2t-6)(x-t) + t^2-6t+7 // 接線の指定通過座標 // x=3, y=-5 // -5=(2t-6)(3-t)+t^2-6t+7 // -5=6t-18-2t^2+6t+t^2-6t+7 // -5=-t^2+6t-11 // t^2-6t+6=0 // t1=(6+√(-6^2-4*6))/2= 4.7320508075688772935274463415059 // t2=(6-√(-6^2-4*6))/2 = 1.2679491924311227064725536584941 //------------------------------------------------------------------------------ // f(t) = (t-3)^2-2 or f(t) = -((t-3)^2-2) // x, y 接線指定通過点座標 // y=(2t-6)(x-t) + t^2-6t+7 // y=2xt-6x-2t^2+6t+t^2-6t+7 // y=-t^2+2xt-6x+7 // t^2-2xt+6x-7+y= 0 or t^2-2xt+6x-7-y= 0 function parabola_tangent(x, y: double; var t1, t2: double): boolean; var a, b, c, d: double; begin if form1.CheckBox1.Checked then y := -y; a := 1; b := - 2 * x; c := 6 * x - 7 + y; d := b * b - 4 * a * c; result := true; if d < 0 then begin result := false; exit; end; d := sqrt(d); t1 := (-b + d) / 2 / a; t2 := (-b - d) / 2 / a; end; //------------------------------------------------------------------------------ // 関数f(x)の計算式 非線形方程式 // f(x) = (x-3)^2 - 2 or f(x) = (x-3)^3 - 2 function f(x: double): double; begin if Form1.Radiobutton2.Checked = true then result := (x - 3) * (x - 3) * (x - 3) - 2 else result := (x - 3) * (x - 3) - 2; if form1.CheckBox1.Checked = true then result := -result; end; // 線分 接線の作図 // x, y 接線通過座標 // a 接線傾斜 procedure line_segment(x, y, a, gl, gr: double); var c : double; begin c := y - a * x; // y = αx + c if gl < x then begin y := a * gl + c; form1.Series2.AddXY(gl, y); end else begin y := a * x + c; form1.Series2.AddXY(x, y); end; y := a * gr + c; Form1.Series2.AddXY(gr, y); Form1.Memo1.Lines.Append('接線式'); if c >= 0 then Form1.Memo1.Lines.Append(' y = ' + floatTostr(a) + 'x + ' + floatTostr(c)) else Form1.Memo1.Lines.Append(' y = ' + floatTostr(a) + 'x ' + floatTostr(c)) end; // 区間[gl,gr]の関数f(x)のグラフを作図 procedure TForm1.chart_graph(gl, gr: double); const n = 100; var x, dx: double; y : double; i : integer; begin series1.clear; // 非線形方程式による曲線 series2.clear; // 接線 series3.clear; // 初期点と検索範囲上限 series4.clear; // 初期間隔と接点、通過点 if gl = gr then exit; // 範囲がゼロだったら終了 dx := (gr - gl) / n; for i := 0 to n do begin x := gl + i * dx; y := f(x); series1.AddXY(x, y); end; end; const orthants = 10; // 二象限跨ぎフラグ // (xp,yp)原点に対する(ex,ey) 角度 - (sx,sy) 角度 左廻り＋ 右廻り－ // sx, sy 非線形線上の前の座標 // ex, ey 非線形線上の後の座標 // xp, yp 接線の通過指定座標 function gradient_difference(sx, sy, ex, ey, xp, yp: double): double; var sf, ef : integer; // 1, 2, 3, 4 象限 xp, yp を原点とするsx,sy ex,eyの象限 sdx, sdy, edx, edy : double; // x, y 差分 salpha, ealpha : double; // 傾斜 begin result := 0; // 通過点(原点)に対する差分 sdx := sx - xp; sdy := sy - yp; edx := ex - xp; edy := ey - yp; sf := 0; ef := 0; // 各座標の原点(xp,yp)に対する象限 if (sdx >= 0) and (sdy >= 0) then sf := 1; if (sdx < 0) and (sdy >= 0) then sf := 2; if (sdx < 0) and (sdy < 0) then sf := 3; if (sdx >= 0) and (sdy < 0) then sf := 4; if (edx >= 0) and (edy >= 0) then ef := 1; if (edx < 0) and (edy >= 0) then ef := 2; if (edx < 0) and (edy < 0) then ef := 3; if (edx >= 0) and (edy < 0) then ef := 4; // 角度rad +y 0～π -y 0～-π salpha := arctan2(sdy, sdx); ealpha := arctan2(edy, edx); // 差分(rad)計算 if sf = 1 then begin if (ef = 1) or (ef = 2) then result := ealpha - salpha; if ef = 4 then result := ealpha - salpha; end; if sf = 2 then begin if (ef = 2) or (ef = 1) then result := ealpha - salpha; if ef = 3 then result := pi + pi + ealpha - salpha; end; if sf = 3 then begin if (ef = 3) or (ef = 4) then result := ealpha - salpha; if ef = 2 then result := ealpha - pi - pi - salpha; end; if sf = 4 then begin if (ef = 4) or (ef = 3) then result := ealpha - salpha; if ef = 1 then result := ealpha - salpha; end; // 差分が二象限分だったら二象限分値 result > pi; if (sf = 1) and (ef = 3) then result := orthants; if (sf = 2) and (ef = 4) then result := orthants; if (sf = 3) and (ef = 1) then result := orthants; if (sf = 4) and (ef = 2) then result := orthants; end; const PLarge = 5500; // 加算ピッチ大フラグ値 XLimit = 5000; // 検索範囲最大値に達したフラグ // 接線の検索と接線の傾き計算 // a 検索開始点 // p 検索ビッチ // xp, yp 接線通過点 // xm 検索終了点 // epsiron 収束判定基準値 // loop 計算数 // ex, ey 接点 // ep 収束判定値 function bisection_method(a, p, xp, yp, xm, epsilon: double; var loop: integer; var ex, ey, ep: double): double; const Maxloop = 1000; var ea : double; sx, sy: double; Fl : boolean; // 移動方向フラグ begin sx := a; // 最初のx座標 sy := f(sx); // 最初のy座標 ex := sx + p; // ピッチ加算次のx座標 ey := f(ex); // 次のy座標 // 座標eとsの xp,yp中心に対する角度差計算 θe-θs ea := gradient_difference(sx, sy, ex, ey, xp, yp); // eaの値が2象限を越えていたら終了 計算ピッチが大きすぎ終了 if ea = orthants then begin loop := PLarge; result := 0; exit; end; Fl := true; // 左廻りフラグ値 if ea < 0 then Fl := false; // 差分ea が－だったら右廻りフラグ値 // 接線端点座標(ex, ey)の近似計算 repeat ex := sx + p; ey := f(ex); // 座標eとsの xp,yp中心に対する角度差計算 θe-θs ea := gradient_difference(sx, sy, ex, ey, xp, yp); if Fl = false then begin // 右廻り if ea > 0 then begin // 左廻りになったら ex := ex - p - p; ey := f(ex); p := p / 4; end; end else begin // 左廻り if ea < 0 then begin // 右廻りになったら ex := ex - p - p; ey := f(ex); p := p / 4; end; end; sx := ex; sy := ey; inc(loop); // ループ数カウント // xの値が検索最大値に到達したか確認 if ex > xm + p then loop := XLimit; ep := epsilon * (abs(ey) + abs(ex) + 1); until (p < ep) or (loop > Maxloop); // 戻り値 接線の傾き計算 result := (ey - yp) / (ex - xp); end; // 計算実行 procedure TForm1.BitBtn1Click(Sender: TObject); var a, p, alpha, xp, yp, gl, gr, xm : double; loop : integer; epsilon, tmp, ex, ey, ep : double; t1, t2 : double; begin a := strTofloat(labelededit3.Text); p := strTofloat(labelededit4.Text); xp := strTofloat(labelededit5.Text); xm := strTofloat(labelededit7.Text); yp := strTofloat(labelededit6.Text); gl := chart1.BottomAxis.Minimum; gr := chart1.BottomAxis.Maximum; memo1.Clear; // 収束判定値計算 epsilon := 1; repeat epsilon := epsilon / 2; tmp := 1 + epsilon; until tmp = 1; epsilon := epsilon * 2; // 非線形方程式の計算 loop := 0; // 近似計算ループ数 ep := epsilon; // 実際の収束判定値 alpha := bisection_method(a, p, xp, yp, xm, epsilon, loop, ex, ey, ep); if loop = XLimit then begin application.MessageBox('初期値と検索範囲の間に接線はありません。','注意',0); memo1.Clear; exit; end; if loop = PLarge then begin application.MessageBox('初期間隔Pを小さくして下さい。','注意',0); memo1.Clear; exit; end; // 接点位置グラフに追加 Series4.AddXY(ex, ey); // 計算結果表示 if Radiobutton2.Checked = true then memo1.Lines.Append('f(x) = (x - 3)^3 - 2') else memo1.Lines.Append('f(x) = (x - 3)^2 - 2'); memo1.Lines.Append('ループ数 ' + intTostr(loop)); memo1.Lines.Append('収束判定値'); memo1.Lines.Append(floatTostr(ep)); memo1.Lines.Append('α= ' + floatTostr(alpha)); memo1.Lines.Append('接点 x = ' + floatTostr(ex)); memo1.Lines.Append('接点 y = ' + floatTostr(ey)); // 接線描画 line_segment(xp, yp, alpha, gl, gr); if Radiobutton2.Checked = false then begin if parabola_tangent(xp, yp, t1, t2) then begin memo1.Lines.Append('t1 = ' + floatTostr(t1)); memo1.Lines.Append('t2 = ' + floatTostr(t2)); end; end; end; // 入力値確認 procedure TForm1.BitBtn2Click(Sender: TObject); var gl, gr, glm, grm : double; ch : integer; a, xp, yp, y, p, xm: double; begin val(labelededit1.Text, glm, ch); if ch <> 0 then begin application.MessageBox('グラフ左最小値に間違いが有ります。','注意',0); exit; end; val(labelededit2.Text, grm, ch); if ch <> 0 then begin application.MessageBox('グラフ右最大値に間違いが有ります。','注意',0); exit; end; val(labelededit3.Text, a, ch); if ch <> 0 then begin application.MessageBox('初期値aに間違いが有ります。','注意',0); exit; end; val(labelededit4.Text, p, ch); if ch <> 0 then begin application.MessageBox('初期間隔pに間違いが有ります。','注意',0); exit; end; if p <=0 then begin application.MessageBox('初期間隔pはゼロ(0)より大きくして下さい。','注意',0); exit; end; val(labelededit7.Text, xm, ch); if ch <> 0 then begin application.MessageBox('検索上限に間違いが有ります。','注意',0); exit; end; if xm - a <= 0 then begin application.MessageBox('検索上限は初期値より大きくして下さい。','注意',0); exit; end; gl := a - (xm - a) / 10; gr := xm + (xm - a) / 10; val(labelededit5.Text, xp, ch); if ch <> 0 then begin application.MessageBox('直線の座標Xに間違いがあります。','注意',0); exit; end; val(labelededit6.Text, yp, ch); if ch <> 0 then begin application.MessageBox('直線の座標Yに間違いがあります。','注意',0); exit; end; if gr < xp then gr := xp + abs(xp / 10); // x方向グラフ上限 if gl > xp then gl := xp - abs(xp / 10); // x方向グラフ下限 if glm < gl then gl := glm; if grm > gr then gr := grm; chart_graph(gl, gr); // 非線形方程式曲線グラフ作図 y := f(a); series3.AddXY(a, y); // 初期位置グラフ追加 y := f(xm); series3.AddXY(xm, y); // 検索上限位置グラフ追加 application.ProcessMessages; if p > (xm - a) / 4 then begin application.MessageBox('初期間隔pが大きすぎます。','注意',0); exit; end; y := f(a + p); series4.AddXY(a + p, y); // 初期計算間隔グラフ追加 series4.AddXY(xp, yp); // 通過指定位置グラフ追加 bitbtn1.Enabled := True; end; procedure TForm1.LabeledEditChange(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.RadioButton1Click(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.RadioButton2Click(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.CheckBox1Click(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.FormCreate(Sender: TObject); var i : integer; begin for i := low(LabeledEdit) to high(LabeledEdit) do begin LabeledEdit[i] := TLabeledEdit(findComponent('LabeledEdit' + intTostr(i))); LabeledEdit[i].OnChange := LabeledEditChange; end; end; end.
Bigdecimalによる計算
演算精度の桁数の半分ぐらいの精度しかでないので、精度が不足する場合は、Bigdecimalでの計算が必要です。
デフォルトの桁数が64桁なので、30桁程度の精度で計算が出来ます。
bigdecimalの組み込みは、ベルヌーイ数その４をbig_mathについては、bigdecimalによるmathを参照してください。
Delphiの多倍長計算にはMPArithもあります、これを使用するのであれば、組み込み方は、ベルヌーイ数その２を参照してください。
MPArithははfunctioｎ形式が使用出来ない為、工夫が必要です。
// 非線形方程式と通過点指定の接線計算 // 曲線上の座標をx方向へ移動し曲線上の点と指定通過点中心角度の回転方向から接線の傾斜の近似値を計算します。 // 計算例は、容易に微分出来る方程式です、あくまでも計算例です。 unit Main; interface uses Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics, Vcl.Controls, Vcl.Forms, Vcl.Dialogs, VclTee.TeeGDIPlus, VCLTee.Series, Vcl.StdCtrls, Vcl.Buttons, Vcl.ExtCtrls, VCLTee.TeEngine, VCLTee.TeeProcs, VCLTee.Chart; type TForm1 = class(TForm) Chart1: TChart; Series1: TLineSeries; LabeledEdit1: TLabeledEdit; LabeledEdit2: TLabeledEdit; LabeledEdit3: TLabeledEdit; LabeledEdit4: TLabeledEdit; LabeledEdit5: TLabeledEdit; LabeledEdit6: TLabeledEdit; LabeledEdit7: TLabeledEdit; BitBtn1: TBitBtn; Memo1: TMemo; Series3: TPointSeries; Series4: TPointSeries; BitBtn2: TBitBtn; Series2: TLineSeries; CheckBox1: TCheckBox; RadioButton1: TRadioButton; RadioButton2: TRadioButton; procedure BitBtn1Click(Sender: TObject); procedure BitBtn2Click(Sender: TObject); procedure CheckBox1Click(Sender: TObject); procedure CheckBox2Click(Sender: TObject); procedure LabeledEditChange(Sender: TObject); procedure FormCreate(Sender: TObject); procedure RadioButton1Click(Sender: TObject); procedure RadioButton2Click(Sender: TObject); private { Private 宣言 } public { Public 宣言 } LabeledEdit : array[1..7] of TLabeledEdit; end; var Form1: TForm1; implementation {$R *.dfm} uses Velthuis.Bigdecimals, Velthuis.bigintegers, big_Math; // f(t) = (t-3)^2-2 or f(t) = -((t-3)^2-2) // x, y 接線指定通過点座標 // y=(2t-6)(x-t) + t^2-6t+7 // y=2xt-6x-2t^2+6t+t^2-6t+7 // y=-t^2+2xt-6x+7 // t^2-2xt+6x-7+y= 0 or t^2-2xt+6x-7-y= 0 function parabola_tangent(x, y: bigdecimal; var t1, t2: bigdecimal): boolean; var a, b, c, d: bigdecimal; one, two, four, six, seven : bigdecimal; begin if Form1.CheckBox1.Checked then y := -y; one := '1'; two := '2'; four := '4'; six := '6'; seven := '7'; a := one; b := - two * x; c := six * x - seven + y; d := b * b - four * a * c; result := true; if d < 0 then begin result := false; exit; end; d := bigdecimal.Sqrt(d, bigdecimal.DefaultPrecision); t1 := (-b + d) / 2 / a; t2 := (-b - d) / 2 / a; end; // 関数f(x)の計算式 非線形方程式 // f(x) = (x-3)^2 - 2 or f(x) = (x-3)^3 - 2 function f(x: bigdecimal): bigdecimal; begin if Form1.Radiobutton2.Checked = true then result := (x - 3) * (x - 3) * (x - 3) - 2 else result := (x - 3) * (x - 3) - 2; if form1.CheckBox1.Checked = true then result := -result; // 乗算により桁数が増加する場合があるのでデフォルト桁数に丸め result := result.RoundToPrecision(bigdecimal.DefaultPrecision); end; // 線分 接線の作図 // x, y 接線通過座標 // a 接線傾斜 procedure line_segment(x, y, a, gl, gr: bigdecimal); var c : bigdecimal; begin c := y - a * x; // y = αx + c if gl < x then begin y := a * gl + c; form1.Series2.AddXY(bigdecimaltodouble(gl), bigdecimaltodouble(y)); end else begin y := a * x + c; form1.Series2.AddXY(bigdecimaltodouble(x), bigdecimaltodouble(y)); end; y := a * gr + c; Form1.Series2.AddXY(bigdecimaltodouble(gr), bigdecimaltodouble(y)); Form1.Memo1.Lines.Append('接線式'); if c >= bigdecimal.Zero then begin Form1.Memo1.Lines.Append(' y = ' + a.ToString + 'x'); Form1.Memo1.Lines.Append(' +' + c.ToString); end else begin Form1.Memo1.Lines.Append(' y = ' + a.ToString + 'x'); Form1.Memo1.Lines.Append(' ' + c.ToString); end; end; // 区間[gl,gr]の関数f(x)のグラフを作図 procedure chart_graph(gl, gr: bigdecimal); const n = 100; var x, dx: bigdecimal; y : double; i : integer; begin Form1.series1.clear; // 非線形方程式による曲線 Form1.series2.clear; // 接線 Form1.series3.clear; // 初期点と検索範囲上限 Form1.series4.clear; // 初期間隔と接点、通過点 if gl = gr then exit; // 範囲がゼロだったら終了 dx := (gr - gl) / n; for i := 0 to n do begin x := gl + i * dx; y := bigdecimaltodouble(f(x)); Form1.series1.AddXY(bigdecimaltodouble(x), y); end; end; const orthants = '10'; // 二象限跨ぎフラグ // (xp,yp)原点に対する(ex,ey) 角度 - (sx,sy) 角度 左廻り＋ 右廻り－ // sx, sy 非線形線上の前の座標 // ex, ey 非線形線上の後の座標 // xp, yp 接線の通過指定座標 function gradient_difference(sx, sy, ex, ey, xp, yp: bigdecimal): bigdecimal; var sf, ef : integer; // 1, 2, 3, 4 象限 xp, yp を原点とするsx,sy ex,eyの象限 sdx, sdy, edx, edy : bigdecimal; // x, y 差分 salpha, ealpha : bigdecimal; // 傾斜 begin result := 0; // 通過点(原点)に対する差分 sdx := sx - xp; sdy := sy - yp; edx := ex - xp; edy := ey - yp; sf := 0; ef := 0; // 各座標の原点(xp,yp)に対する象限 if (sdx >= 0) and (sdy >= 0) then sf := 1; if (sdx < 0) and (sdy >= 0) then sf := 2; if (sdx < 0) and (sdy < 0) then sf := 3; if (sdx >= 0) and (sdy < 0) then sf := 4; if (edx >= 0) and (edy >= 0) then ef := 1; if (edx < 0) and (edy >= 0) then ef := 2; if (edx < 0) and (edy < 0) then ef := 3; if (edx >= 0) and (edy < 0) then ef := 4; // 角度rad +y 0～π -y 0～-π salpha := arctan2_big(sdy, sdx); ealpha := arctan2_big(edy, edx); // 差分(rad)計算 if sf = 1 then begin if (ef = 1) or (ef = 2) then result := ealpha - salpha; if ef = 4 then result := ealpha - salpha; end; if sf = 2 then begin if (ef = 2) or (ef = 1) then result := ealpha - salpha; if ef = 3 then result := pi_big + pi_big + ealpha - salpha; end; if sf = 3 then begin if (ef = 3) or (ef = 4) then result := ealpha - salpha; if ef = 2 then result := ealpha - pi_big - pi_big - salpha; end; if sf = 4 then begin if (ef = 4) or (ef = 3) then result := ealpha - salpha; if ef = 1 then result := ealpha - salpha; end; // 差分が二象限分だったら二象限分値 result > pi; if (sf = 1) and (ef = 3) then result := orthants; if (sf = 2) and (ef = 4) then result := orthants; if (sf = 3) and (ef = 1) then result := orthants; if (sf = 4) and (ef = 2) then result := orthants; end; const PLarge = 5500; // 加算ピッチ大フラグ値 XLimit = 6000; // 検索範囲最大値に達したフラグ // 接線の検索と接線の傾き計算 // a 検索開始点 // p 検索ビッチ // xp, yp 接線通過点 // xm 検索終了点 // epsiron 収束判定基準値 // loop 計算数 // ex, ey 接点 // ep 収束判定値 function bisection_method(a, p, xp, yp, xm, epsilon: bigdecimal; var loop: integer; var ex, ey, ep: bigdecimal): bigdecimal; const Maxloop = 5000; var ea : bigdecimal; sx, sy, four: bigdecimal; defa, defb, orth : bigdecimal; Fl : boolean; // 移動方向フラグ begin orth := orthants; four := '4'; sx := a; // 最初のx座標 sy := f(sx); // 最初のy座標 ex := sx + p; // ピッチ加算次のx座標 ey := f(ex); // 次のy座標 // 座標eとsの xp,yp中心に対する角度差計算 θe-θs ea := gradient_difference(sx, sy, ex, ey, xp, yp); // eaの値が2象限を越えていたら終了 計算ピッチが大きすぎ終了 if ea = orth then begin loop := PLarge; result := 0; exit; end; Fl := true; // 左廻りフラグ値 if ea < 0 then Fl := false; // 差分ea が－だったら右廻りフラグ値 // 接線端点座標(ex, ey)の近似計算 repeat ex := sx + p; ey := f(ex); // 座標eとsの xp,yp中心に対する角度差計算 θe-θs ea := gradient_difference(sx, sy, ex, ey, xp, yp); if Fl = false then begin // 右廻り if ea > 0 then begin // 左廻りになったら ex := ex - p - p; ey := f(ex); p := p / four; end; end else begin // 左廻り if ea < 0 then begin // 右廻りになったら ex := ex - p - p; ey := f(ex); p := p / four; end; end; sx := ex; sy := ey; inc(loop); // ループ数カウント // xの値が検索最大値に到達したか確認 if ex > xm + p then loop := XLimit; ep := epsilon * (ey.Abs(ey) + ex.abs(ex) + 1); until (p < ep) or (loop > Maxloop); // 戻り値 接線の傾き計算 defa := ey - yp; defb := ex - xp; // defa := defa.RoundToPrecision(bigdecimal.DefaultPrecision); // defb := defb.RoundToPrecision(bigdecimal.DefaultPrecision); result := defa / defb; // result := bigdecimal.Divide(defa, defb); end; // 計算実行 procedure TForm1.BitBtn1Click(Sender: TObject); const DIG = 40; // Bigdecimal 表示桁数 var a, p, alpha, xp, yp, gr, gl, xm: bigdecimal; loop, prec : integer; epsilon, ex, ey, ep : bigdecimal; estr : string; t1, t2 : bigdecimal; begin a := labelededit3.Text; p := labelededit4.Text; xp := labelededit5.Text; xm := labelededit7.Text; yp := labelededit6.Text; gl := chart1.BottomAxis.Minimum; gr := chart1.BottomAxis.Maximum; memo1.Clear; // 収束判定値計算 prec := bigdecimal.DefaultPrecision; prec := prec - 1; estr := '1e-' + intTostr(prec); epsilon := estr; // 非線形方程式の計算 loop := 0; // 近似計算ループ数 ep := epsilon; // 実際の収束判定値 alpha := bisection_method(a, p, xp, yp, xm, epsilon, loop, ex, ey, ep); if loop = XLimit then begin application.MessageBox('初期値と検索範囲の間に接線はありません。','注意',0); memo1.Clear; exit; end; if loop = PLarge then begin application.MessageBox('初期間隔Pを小さくして下さい。','注意',0); memo1.Clear; exit; end; // 接点位置グラフに追加 Series4.AddXY(bigdecimaltodouble(ex), bigdecimaltodouble(ey)); // 計算結果表示 if Radiobutton2.Checked = true then memo1.Lines.Append('f(x) = (x - 3)^3 - 2') else memo1.Lines.Append('f(x) = (x - 3)^2 - 2'); memo1.Lines.Append('ループ数 ' + intTostr(loop)); memo1.Lines.Append('収束判定値'); ep := ep.RoundToPrecision(DIG); memo1.Lines.Append(ep.ToString); alpha := alpha.RoundToPrecision(DIG); memo1.Lines.Append('α= ' + alpha.ToString); ex := ex.RoundToPrecision(DIG); ey := ey.RoundToPrecision(DIG); memo1.Lines.Append('接点 x = ' + ex.ToString); memo1.Lines.Append('接点 y = ' + ey.ToString); // 接線描画 line_segment(xp, yp, alpha, gl, gr); if Radiobutton2.Checked = false then begin if parabola_tangent(xp, yp, t1, t2) then begin memo1.Lines.Append('t1 = ' + t1.ToString); memo1.Lines.Append('t2 = ' + t2.ToString); end; end; end; // 入力値確認 procedure TForm1.BitBtn2Click(Sender: TObject); var gl, gr, glm, grm : double; ch : integer; a, xp, yp, y, p, xm: double; begin val(labelededit1.Text, glm, ch); if ch <> 0 then begin application.MessageBox('グラフ左最小値に間違いが有ります。','注意',0); exit; end; val(labelededit2.Text, grm, ch); if ch <> 0 then begin application.MessageBox('グラフ右最大値に間違いが有ります。','注意',0); exit; end; val(labelededit3.Text, a, ch); if ch <> 0 then begin application.MessageBox('初期値aに間違いが有ります。','注意',0); exit; end; val(labelededit4.Text, p, ch); if ch <> 0 then begin application.MessageBox('初期間隔pに間違いが有ります。','注意',0); exit; end; if p <=0 then begin application.MessageBox('初期間隔pはゼロ(0)より大きくして下さい。','注意',0); exit; end; val(labelededit7.Text, xm, ch); if ch <> 0 then begin application.MessageBox('検索上限に間違いが有ります。','注意',0); exit; end; if a >= xm then begin application.MessageBox('検索上限は初期値aより大きくしてください。','注意',0); exit; end; val(labelededit5.Text, xp, ch); if ch <> 0 then begin application.MessageBox('直線の座標Xに間違いがあります。','注意',0); exit; end; val(labelededit6.Text, yp, ch); if ch <> 0 then begin application.MessageBox('直線の座標Yに間違いがあります。','注意',0); exit; end; gl := a - (xm - a) / 10; gr := xm + (xm - a) / 10; if xp < gl then gl := xp - abs(xp / 10); if xp > gr then gr := xp + abs(xp / 10); if glm < gl then gl := glm; if grm > gr then gr := grm; chart_graph(gl, gr); // 非線形方程式曲線グラフ作図 y := bigdecimaltodouble(f(a)); series3.AddXY(a, y); // 初期位置グラフ追加 y := bigdecimaltodouble(f(xm)); series3.AddXY(xm, y); // 検索上限位置グラフ追加 application.ProcessMessages; if p > (xm - a) / 4 then begin application.MessageBox('初期間隔pが大きすぎます。','注意',0); exit; end; y := bigdecimaltodouble(f(a + p)); series4.AddXY(a + p, y); // 初期計算間隔グラフ追加 series4.AddXY(xp, yp); // 通過指定位置グラフ追加 bitbtn1.Enabled := True; end; procedure TForm1.LabeledEditChange(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.RadioButton1Click(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.RadioButton2Click(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.CheckBox1Click(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.CheckBox2Click(Sender: TObject); begin bitbtn1.Enabled := false; end; procedure TForm1.FormCreate(Sender: TObject); var i : integer; begin for i := low(LabeledEdit) to high(LabeledEdit) do begin LabeledEdit[i] := TLabeledEdit(findComponent('LabeledEdit' + intTostr(i))); LabeledEdit[i].OnChange := LabeledEditChange; end; end; end.
tangent_line_point_double.zip.zip
三角関数、逆三角関数、その他関数、 連立一次方程式の解法 に戻る