Skip to content

Commit

Permalink
#65 update solve operator
Browse files Browse the repository at this point in the history
  • Loading branch information
shapoco committed Jun 2, 2023
1 parent bd7aed8 commit e348d2e
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 28 deletions.
78 changes: 57 additions & 21 deletions Calctus/Model/Expressions/SolveExpr.cs
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,14 @@
using Shapoco.Calctus.Model.Mathematics;
using Shapoco.Calctus.Model.Parsers;
using Shapoco.Calctus.Model.Types;
using Shapoco.Calctus.Model.Formats;

namespace Shapoco.Calctus.Model.Expressions {
class SolveExpr : Expr {
public readonly Expr Equation;
public readonly Token Variant;
public readonly Expr ParamMin;
public readonly Expr ParamMax;
public readonly Expr Param0;
public readonly Expr Param1;

public SolveExpr(Token keyword, Expr equation, Token variant, Expr paramMin, Expr paramMax) : base(keyword) {
// 等式が与えられた場合は減算にすり替える
Expand All @@ -23,33 +24,68 @@ public SolveExpr(Token keyword, Expr equation, Token variant, Expr paramMin, Exp

Equation = equation;
Variant = variant;
ParamMin = paramMin;
ParamMax = paramMax;
Param0 = paramMin;
Param1 = paramMax;
}

protected override Val OnEval(EvalContext e) {
// 初期値の範囲
var pMinVal = ParamMin != null ? ParamMin.Eval(e) : new RealVal(-10m);
var pMaxVal = ParamMax != null ? ParamMax.Eval(e) : new RealVal(10m);
var pMin = (decimal)pMinVal.AsReal;
var pMax = (decimal)pMaxVal.AsReal;
if (pMin >= pMax) {
throw new ArgumentException("Invalid parameter range.");
}
// パラメータの評価
var paramVal0 = Param0 != null ? Param0.Eval(e) : null;
var paramVal1 = Param1 != null ? Param1.Eval(e) : null;
var formatHint = new FormatHint(NumberFormatter.CStyleReal);

// 数値微分用の定数
var h = Math.Max(1e-20m, (pMax - pMin) / 1e6m);
var initVals = new List<decimal>();
var h = 1e-5m;
var tol = 1e-10m;
var pMin = decimal.MinValue;
var pMax = decimal.MaxValue;
const int N = 100;

// 収束条件
var tol = Math.Max(1e-25m, (pMax - pMin) / 1e12m);
if (paramVal0 != null && paramVal1 != null) {
// 初期値の範囲が指定された場合
// --> 範囲を N 分割して初期値を決定
pMin = (decimal)paramVal0.AsReal;
pMax = (decimal)paramVal1.AsReal;
if (pMin >= pMax) {
throw new ArgumentException("Invalid parameter range.");
}
for (int i = 0; i <= N; i++) {
initVals.Add(pMin + (pMax - pMin) * i / N);
}
h = Math.Max(1e-16m, (pMax - pMin) / 1e6m);
tol = Math.Max(1e-20m, (pMax - pMin) / 1e12m);
formatHint = paramVal0.FormatHint;
}
else if (paramVal0 != null) {
// 初期値が直接指定された場合
if (paramVal0 is ArrayVal) {
// 初期値が配列で指定された場合
foreach(var val in paramVal0.AsRealArray) {
initVals.Add(val);
}
var paramValArray = (Val[])paramVal0.Raw;
if (paramValArray.Length > 0) {
formatHint = paramValArray[0].FormatHint;
}
}
else {
// 初期値がスカラ値で指定された場合
initVals.Add(paramVal0.AsReal);
formatHint = paramVal0.FormatHint;
}
}
else {
// 引数が指定されなかった場合
for (int i = 0; i <= N; i++) {
initVals.Add(-10 + 20 * i / N);
}
}

// ニュートン法
var scope = new EvalContext(e);
var results = new List<decimal>();
const int N = 100;
for (int i = 0; i < N; i++) {
foreach (var init in initVals) {
try {
var init = pMin + (pMax - pMin) * i / N;
if (newtonMethod(scope, init, pMin, pMax, h, tol, out decimal r)) {
if (!results.Any(p => Math.Abs(p - r) < tol * 2)) {
results.Add(r);
Expand All @@ -60,11 +96,11 @@ protected override Val OnEval(EvalContext e) {
}

if (results.Count == 1) {
return new RealVal(results[0]).Format(pMinVal.FormatHint);
return new RealVal(results[0]).Format(formatHint);
}
else {
results.Sort();
return new ArrayVal(results.Select(p => new RealVal(p).Format(pMinVal.FormatHint)).ToArray());
return new ArrayVal(results.Select(p => new RealVal(p).Format(formatHint)).ToArray());
}
}

Expand Down
13 changes: 7 additions & 6 deletions Calctus/Model/Parsers/Parser.cs
Original file line number Diff line number Diff line change
Expand Up @@ -174,15 +174,16 @@ public Expr Solve(Token first) {
var equation = Expr(false);
Expect(",");
Expect(TokenType.Word, out Token variant);
Expr min = null;
Expr max = null;
Expr param0 = null;
Expr param1 = null;
if (ReadIf(",")) {
min = Expr(false);
Expect(",");
max = Expr(false);
param0 = Expr(false);
if (ReadIf(",")) {
param1 = Expr(false);
}
}
Expect(")");
return new SolveExpr(first, equation, variant, min, max);
return new SolveExpr(first, equation, variant, param0, param1);
}

public Token Peek() {
Expand Down
2 changes: 1 addition & 1 deletion Calctus/UI/Sheets/SheetView.cs
Original file line number Diff line number Diff line change
Expand Up @@ -819,7 +819,7 @@ private void processRecalcRequest() {
list.Add(new InputCandidate(BoolVal.TrueKeyword, BoolVal.TrueKeyword, "true value", false));
list.Add(new InputCandidate(BoolVal.FalseKeyword, BoolVal.FalseKeyword, "false value", false));
list.Add(new InputCandidate("def", "def", "user function definition", false));
list.Add(new InputCandidate("solve", "solve(expr,var,min=-10,max=10)", "Solves the equation using Newton's method.", true));
list.Add(new InputCandidate("solve", "solve(expr,var,a,b)", "Solves the equation using Newton's method.", true));
_inputCandidates = list.OrderBy(p => p.Id).ToArray();
_recalcRequested = false;
}
Expand Down

0 comments on commit e348d2e

Please sign in to comment.