作者Lucemia (生の直感、死の予感)
看板C_Sharp
標題[心得] C# 代數處理的類
時間Thu Jul 16 14:35:17 2009
今天臨時需要用到解方程式，牛頓法
找不到C#的類別。
python的 sympy又跑太慢了,
只好自己硬幹了一個能跑牛頓法的。
請問C#有這類的lib嗎?
// 用法
// 定義變數
Furmula x = new Furmula("x");
Furmula y = new Furmula("y");
Furmula z = new Furmula("z");
// 定義方程式
Furmula f = x*x + y*y+z*z + x*y*z;
// 變數代入
double v =f.subs(x,3).subs(y,2).subs(z,1).Value;
// 微分
Furmula df = f.Diff(x);
// -- source code --
class Furmula
{
// Furmula
public String mName = "";
public bool mIsValue = false;
public double mValue = 0;
//
public OPERATORS mOperator = OPERATORS.NULL;
public Furmula mLeftChild = null;
public Furmula mRightChild = null;
public enum OPERATORS { ADD, SUB, MUL, DIV, SIN, COS, TAN, POW, NULL };
public Furmula(String name)
{
mName = name;
}
public Furmula(double value)
{
mValue = value;
mIsValue = true;
}
public Furmula(Furmula b)
{
mName = b.mName;
mIsValue = b.mIsValue;
mValue = b.mValue;
mOperator = b.mOperator;
mLeftChild = b.mLeftChild;
mRightChild = b.mRightChild;
}
public Furmula(OPERATORS oper, Furmula LeftChild, Furmula RightChild)
{
mOperator = oper;
mLeftChild = LeftChild;
mRightChild = RightChild;
}
public static Furmula operator +(Furmula a, Furmula b)
{
if (a.mIsValue && b.mIsValue)
return a.mValue + b.mValue;
else if (a == 0)
return new Furmula(b);
else if (b == 0)
return new Furmula(a);
return new Furmula(OPERATORS.ADD, a, b);
}
public static Furmula operator -(Furmula a, Furmula b)
{
if (a.mIsValue && b.mIsValue)
return a.mValue - b.mValue;
else if (a == 0)
return -1 * new Furmula(b);
else if (b == 0)
return new Furmula(a);
return new Furmula(OPERATORS.SUB, a, b);
}
public static Furmula operator *(Furmula a, Furmula b)
{
if (a.mIsValue && b.mIsValue)
return a.mValue * b.mValue;
else if (a == 0)
return 0;
else if (b == 0)
return 0;
else if (a == 1)
return new Furmula(b);
else if (b == 1)
return new Furmula(a);
if (a.mIsValue && b.mIsValue)
return a.mValue * b.mValue;
return new Furmula(OPERATORS.MUL, a, b);
}
public static Furmula operator /(Furmula a, Furmula b)
{
if (a.mIsValue && b.mIsValue)
return a.mValue / b.mValue;
else if (b == 0)
throw new Exception();
else if (a == 0)
return 0;
else if (b == 1)
return new Furmula(a);
return new Furmula(OPERATORS.DIV, a, b);
}
public static Furmula Pow(Furmula a, Furmula b)
{
if (a.mIsValue && b.mIsValue)
return Math.Pow(a.mValue, b.mValue);
else if (b == 0)
return 1;
else if (a == 0)
return 0;
else if (a == 1)
return 1;
return new Furmula(OPERATORS.POW, a, b);
}
public static Furmula Sin(Furmula a)
{
if (a.mIsValue)
return Math.Sin(a.mValue);
return new Furmula(OPERATORS.SIN, a, null);
}
public static Furmula Cos(Furmula a)
{
if (a.mIsValue)
return Math.Cos(a.mValue);
return new Furmula(OPERATORS.COS, a, null);
}
public static Furmula Tan(Furmula a)
{
if (a.mIsValue)
return Math.Tan(a.mValue);
return new Furmula(OPERATORS.TAN, a, null);
}
public static bool operator ==(Furmula a, double b)
{
if (a.mIsValue && a.mValue == b)
return true;
return false;
}
public static bool operator !=(Furmula a, double b)
{
return !(a == b);
}
public static implicit operator Furmula(double v)
{
return new Furmula(v);
}
public double Value
{
get
{
if (mIsValue)
{
return mValue;
}
throw new Exception();
}
}
public Furmula Subs(Dictionary<Furmula, double> values)
{
if (mOperator == OPERATORS.NULL)
{
foreach (var v in values)
{
if (mName == v.Key.mName)
{
return new Furmula(v.Value);
}
}
return new Furmula(this);
}
else
{
Furmula NewLeft = null, NewRight = null;
if (mLeftChild != null)
{
NewLeft = mLeftChild.Subs(values);
}
if (mRightChild != null)
{
NewRight= mRightChild.Subs(values);
}
switch (mOperator)
{
case OPERATORS.ADD:
return NewLeft + NewRight;
case OPERATORS.SUB:
return NewLeft - NewRight;
case OPERATORS.MUL:
return NewLeft * NewRight;
case OPERATORS.DIV:
return NewLeft / NewRight;
case OPERATORS.SIN:
return Sin(NewLeft);
case OPERATORS.COS:
return Cos(NewLeft);
case OPERATORS.TAN:
return Tan(NewLeft);
case OPERATORS.POW:
return Pow(NewLeft, NewRight);
}
throw new Exception();
}
}
public Furmula Subs(String name, double value)
{
if (mOperator == OPERATORS.NULL)
{
if (mName == name)
{
return new Furmula(value);
}
else
{
return new Furmula(this);
}
}
else
{
Furmula NewLeft = null, NewRight = null;
if (mLeftChild != null)
{
NewLeft= mLeftChild.Subs(name, value);
}
if (mRightChild != null)
{
NewRight= mRightChild.Subs(name, value);
}
switch (mOperator)
{
case OPERATORS.ADD:
return NewLeft + NewRight;
case OPERATORS.SUB:
return NewLeft - NewRight;
case OPERATORS.MUL:
return NewLeft * NewRight;
case OPERATORS.DIV:
return NewLeft / NewRight;
case OPERATORS.SIN:
return Sin(NewLeft);
case OPERATORS.COS:
return Cos(NewLeft);
case OPERATORS.TAN:
return Tan(NewLeft);
case OPERATORS.POW:
return Pow(NewLeft, NewRight);
}
throw new Exception();
}
}
public Furmula Diff(Furmula furmula)
{
return Diff(furmula.mName);
}
public Furmula Diff(String name)
{
if (mOperator == OPERATORS.NULL)
{
if (mName == name)
{
return new Furmula(1);
}
else
{
return new Furmula(0);
}
}
else
{
switch (mOperator)
{
case OPERATORS.ADD:
return mLeftChild.Diff(name) + mRightChild.Diff(name);
case OPERATORS.SUB:
return mLeftChild.Diff(name) - mRightChild.Diff(name);
case OPERATORS.MUL:
return mLeftChild.Diff(name) * mRightChild + mLeftChild * mRightChild.Diff(name);
case OPERATORS.DIV:
return ((mLeftChild.Diff(name) * mRightChild) - (mLeftChild * mRightChild.Diff(name))) / Pow(mRightChild, 2);
case OPERATORS.SIN:
return Cos(mLeftChild) * mLeftChild.Diff(name);
case OPERATORS.COS:
return -1 * Sin(mLeftChild) * mLeftChild.Diff(name);
case OPERATORS.TAN:
return mLeftChild.Diff(name) / Pow(Cos(mLeftChild), 2);
case OPERATORS.POW:
return mRightChild * Pow(mLeftChild, mRightChild - 1) * mLeftChild.Diff(name);
}
}
throw new Exception();
}
public override string ToString()
{
switch(mOperator){
case OPERATORS.NULL:
if (mIsValue)
return mValue.ToString();
else
return mName;
case OPERATORS.ADD:
return String.Format("({0}+{1})", mLeftChild, mRightChild);
case OPERATORS.SUB:
return String.Format("({0}-{1})", mLeftChild, mRightChild);
case OPERATORS.MUL:
return String.Format("({0}*{1})", mLeftChild, mRightChild);
case OPERATORS.DIV:
return String.Format("({0}/{1})", mLeftChild, mRightChild);
case OPERATORS.POW:
return String.Format("({0}**{1})", mLeftChild, mRightChild);
case OPERATORS.SIN:
return String.Format("Sin({0})", mLeftChild);
case OPERATORS.COS:
return String.Format("Cos({0})", mLeftChild);
case OPERATORS.TAN:
return String.Format("Tan({0})", mLeftChild);
}
throw new Exception();
}
public static Dictionary<Furmula, double> Newton(Furmula f, Dictionary<Furmula, double> Guass)
{
List<Furmula> vars = new List<Furmula>(Guass.Keys);
int K = vars.Count;
Furmula[] Gradient = new Furmula[K];
Furmula[,] Hessian = new Furmula[K, K];
for (int i = 0; i < K; i++)
Gradient[i] = f.Diff(vars[i]);
for (int i = 0; i < K; i++)
{
for (int j = 0; j < K; j++)
{
Hessian[i, j] = f.Diff(vars[i]).Diff(vars[j]);
}
}
CSML.Matrix Xn = new CSML.Matrix(K, 1);
CSML.Matrix Fn = new CSML.Matrix(K, 1);
CSML.Matrix Hn = new CSML.Matrix(K, K);
System.IO.StreamWriter writer = new System.IO.StreamWriter("iter.log");
for (int iter = 0; iter < 100; iter++)
{
writer.WriteLine("{0} {1}", iter, f.Subs(Guass));
writer.Flush();
for (int i = 0; i < K; i++)
{
Xn[i + 1, 1] = new Complex(Guass[vars[i]]);
}
for (int i = 0; i < K; i++)
{
Fn[i + 1, 1] = new Complex(Gradient[i].Subs(Guass).Value);
}
for (int i = 0; i < K; i++)
{
for (int j = 0; j < K; j++)
{
Hn[i + 1, j + 1] = new Complex(Hessian[i, j].Subs(Guass).Value);
}
}
Xn = Xn - Hn.Inverse() * Fn;
for (int i = 0; i < K; i++)
{
Guass[vars[i]] = Xn[i + 1, 1].Re;
}
}
return Guass;
}
}
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 154.20.172.64
1F：推 Steven0422:IMSL for .Net 07/16 16:34
2F：推 pyrochlore:Open Source的有Math.Net 07/17 14:08
3F：→ Lucemia:感謝 07/18 00:32