385. Dual numbers

נחזור קצת לטיפים בנושאים מתמטיים:

אחד הדברים שלומדים בתיכון זה כיצד לגזור פונקציה. לעתים זה יכול להיות שימושי לגזור פונקציה בקוד.

למשל, אלגוריתם ניוטון רפסון מאפשר לנו למצוא אפסים של פונקציה בצורה יחסית מהירה. (דיברנו על זה קצת בטיפ 341).

האלגוריתם הזה דורש מאיתנו לדעת בכל שלב מהי הנגזרת של הפונקציה בנקודה הנוכחית.

---

קיימים מספר סוגי גזירות של פונקציה:

• גזירה נומרית – מאפשרת לנו לקבל קירוב לערך של הנגזרת של הפונקציה בנקודה מסוימת
• גזירה סימבולית – מאפשרת לנו להתייחס לפונקציה כביטוי מתמטי, וע”י פעולת הגזירה לקבל ביטוי מתמטי אחר.

על גזירה סימבולית לא נדבר הפעם, אך רק נציין שיש לגזירה זו קשר לExpression Trees שראינו ב.net 3.5 (טיפים מספר 173-181)

כעת גם גזירה נומרית מתחלקת למספר סוגים:

הסוג הראשון הוא סוג שנותן קירוב לנגזרת של הפונקציה בנקודה ע”י כל מיני חישובים מתמטיים המסתמכים על אנליזה מתמטית ונותן גם הערכה לשגיאה. נוסחאות כאלה חוקרים בקורסים במתמטיקה בנושאי אנליזה נומרית.

הסוג השני הוא סוג שנותן את הערך הנגזרת של הפונקציה בצורה מדויקת – עד כדי חישוב: למה הכוונה? הכוונה היא שאם נסתכל על האלגוריתם המתמטי שבו מחושבת הנוסחה, נראה שהנוסחה נכונה לחלוטין, אבל למחשב שבו אנחנו נריץ את האלגוריתם יש דיוק מסוים בפונקציות בהן הוא משתמש, ולכן החישוב מושפע ממנו.

נדבר על הסוג השני כעת.

---

אז איך אפשר לחשב נגזרת של פונקציה בצורה מדויקת? שיטה אחת נקראת Dual Numbers והיא מוגדרת ככה:

נגדיר מספר דואלי בתור זוג סדור $ (x, y)$ של מספרים ממשיים. נגדיר עליו את הפעולות הבאות:

• $ \left(u, u’ \right) + \left( v, v’ \right) = \left(u + v, u’ + v’ \right) $
• $ \left(u, u’ \right) - \left( v, v’ \right) = \left(u - v, u’ - v’ \right) $
• $ \left(u, u’ \right) \cdot \left( v, v’ \right) = \left(u \cdot v, u \cdot v’ + v \cdot u’ \right) $
• $ \displaystyle{\frac { \left(u, u’ \right) }{ \left( v, v’ \right) } = \left( \frac{u}{v}, \frac{u’ \cdot v - u \cdot v’}{v^2} \right) } $

בנוסף נגדיר לפונקציה גזירה f פונקציה דואלית F $ F\left(u, u’ \right) = \left(f(u), f’(u) \cdot u’ \right) $, למשל

• $ \sin \left(u, u’ \right) = \left(\sin u, \cos u \cdot u’ \right) $
• $ \cos\left(u, u’ \right) = \left(\cos u, -\sin u \cdot u’ \right) $
• $ \exp\left(u, u’ \right) = \left(\exp u, \exp u \cdot u’ \right) $

בנוסף נתייחס לקבוע $ C $ בתור הערך $ \left( C, 0 \right) $

---

למה זה טוב? נגדיר ואנחנו רוצים לחשב את הנגזרת של הפונקציה $ \displaystyle{\frac{\sin x}{x^2 + 1}} $:

נתחיל עם הזוג $ \left( x, 1 \right) $ (כי הנגזרת של הפונקציה x בכל נקודה היא 1).

נבצע עליה פעולת $ \sin $ ונקבל $ \left( \sin x, \cos x \right) $.

בנוסף נבצע $ \left( x, 1 \right) \cdot \left( x, 1 \right) = \left( x \cdot x, x \cdot 1 + 1 \cdot x \right) = \left( x^2, 2x \right) $ ונוסיף לזה $ \left(1, 0\right) $ ונקבל $ \left(x^2 + 1, 2x \right) $. כעת נבצע את הפעולה $ \frac{\left(\sin x, \cos x \right)}{\left(x^2 + 1, 2x \right)}$ ונקבל $ \displaystyle{ \left(\frac{\sin x}{x^2 + 1}, \frac{\cos x \cdot (x^2 + 1) - \sin x \cdot (2x)}{(x^2 + 1)^2} \right) } $

שימו לב שקיבלנו שהנגזרת של $ \displaystyle{\frac{\sin x}{x^2 + 1}} $ היא $ \displaystyle{ \frac{\cos x \cdot (x^2 + 1) - \sin x \cdot (2x)}{(x^2 + 1)^2} } $

קיבלנו את זה בלי להשתמש ישירות בנוסחאות של התיכון, אלא רק בנוסחאות שרשמתי מעלה.

בעצם הנוסחאות שרשמתי מעלה הן נוסחאות שמתייחסות לכל מספר דואלי בתור זוג סדור: ערך של פונקציה בנקודה וערך של הנגזרת של הפונקציה בנקודה.

הפעולות הבינאריות והאונאריות שהוגדרו למעלה, פועלות “כרגיל” על החלק שמייצג את הערך של הפונקציה בנקודה, ומכבדות את חוקי הגזירה בחלק שמייצג את הנגזרת בנקודה.

---

אז איך נשתמש בזה בתכנות? קל לכתוב מחלקה שתממש את הפונקציונאליות, למשל:

public class DualNumber
{
    #region Data Members

    private double m_FunctionValue;
    private double m_DerivativeValue;

    private static readonly DualNumber m_Zero = new DualNumber(0, 0);
    private static readonly DualNumber m_One = new DualNumber(1, 0);

    #endregion

    #region Constructors

    public DualNumber(double functionValue,
                      double derivativeValue)
    {
        m_FunctionValue = functionValue;
        m_DerivativeValue = derivativeValue;
    }

    public DualNumber(DualNumber source) :
        this(source.FunctionValue,
             source.DerivativeValue)
    {
    }

    #endregion

    #region Constants
    public static DualNumber ONE
    {
        get
        {
            return m_One;
        }
    }

    public static DualNumber ZERO
    {
        get
        {
            return m_Zero;
        }
    }

    #endregion
    
    #region Binary Methods
    
    private DualNumber Add(DualNumber givenDual)
    {
        return new DualNumber(this.FunctionValue + givenDual.FunctionValue,
                              this.DerivativeValue + givenDual.DerivativeValue);
    }

    private DualNumber Sub(DualNumber givenDual)
    {
        return new DualNumber(this.FunctionValue - givenDual.FunctionValue,
                              this.DerivativeValue - givenDual.DerivativeValue);
    }

    private DualNumber Mul(DualNumber givenDual)
    {
        return new DualNumber(this.FunctionValue * givenDual.FunctionValue,
                              this.DerivativeValue * givenDual.FunctionValue+
                              this.FunctionValue * givenDual.DerivativeValue);
    }

    private DualNumber Div(DualNumber givenDual)
    {
        return new DualNumber(this.FunctionValue / givenDual.FunctionValue,
                             (this.DerivativeValue * givenDual.FunctionValue -
                              this.FunctionValue * givenDual.DerivativeValue)/
                              (givenDual.FunctionValue * givenDual.FunctionValue));
    }
    
    private DualNumber Power(int exp)
    {
        if (exp >= 0)
        {
            return this.PositivePower(exp);
        }
        else
        {
            return (this.Opposite.PositivePower(-exp));
        }
    }

    private DualNumber Opposite
    {
        get
        {
            return new DualNumber(1/this.FunctionValue,
                                  -this.DerivativeValue/(this.FunctionValue * this.FunctionValue));
        }
    }

    private DualNumber PositivePower(int exp)
    {
        DualNumber powerResult = DualNumber.ONE;
        DualNumber loopPower = new DualNumber(this);
        int loopCurrentPower = exp;
        while (loopCurrentPower != 0)
        {
            if (loopCurrentPower % 2 == 1)
            {
                powerResult = powerResult.Mul(loopPower);
            }

            loopPower = loopPower.Mul(loopPower);
            loopCurrentPower /= 2;
        }

        return powerResult;
    }

    public override bool Equals(object obj)
    {
        if (!(obj is DualNumber))
        {
            return false;
        }
        else
        {
            DualNumber comparedDual = (DualNumber)obj;

            return (this.FunctionValue.Equals(comparedDual.FunctionValue) &&
                    this.DerivativeValue.Equals(comparedDual.DerivativeValue));
        }
    }

    public override int GetHashCode()
    {
        return (this.FunctionValue.GetHashCode() ^
                this.DerivativeValue.GetHashCode());
    }

    #endregion

    #region Unary Methods

    public double FunctionValue
    {
        get
        {
            return m_FunctionValue;
        }
    }
    public double DerivativeValue
    {
        get
        {
            return m_DerivativeValue;
        }
    }

    #endregion

    #region Operators

    public static DualNumber operator +(DualNumber z1, DualNumber z2)
    {
        return z1.Add(z2);
    }

    public static DualNumber operator -(DualNumber z1, DualNumber z2)
    {
        return z1.Sub(z2);
    }

    public static DualNumber operator *(DualNumber z1, DualNumber z2)
    {
        return z1.Mul(z2);
    }

    public static DualNumber operator /(DualNumber z1, DualNumber z2)
    {
        return z1.Div(z2);
    }

    public static DualNumber operator ^(DualNumber z, int exp)
    {
        return z.Power(exp);
    }

    public static bool operator ==(DualNumber z1, DualNumber z2)
    {
        return z1.Equals(z2);
    }

    public static bool operator !=(DualNumber z1, DualNumber z2)
    {
        return (!z1.Equals(z2));
    }

    #endregion
    
    #region Elementary functions
    
    public static DualNumber Cos(DualNumber argument)
    {
        return new DualNumber(Math.Cos(argument.FunctionValue),
                              -Math.Sin(argument.FunctionValue)*
                              argument.DerivativeValue);
    }

    public static DualNumber Sin(DualNumber argument)
    {
        return new DualNumber(Math.Sin(argument.FunctionValue),
                              Math.Cos(argument.FunctionValue) *
                              argument.DerivativeValue);
    }

    public static DualNumber Exp(DualNumber argument)
    {
        return new DualNumber(Math.Exp(argument.FunctionValue),
                              Math.Exp(argument.FunctionValue) *
                              argument.DerivativeValue);
    }

    #endregion
}

כעת נוכל לכתוב פונקציה שמחשבת את הנגזרת של הפונקציה $ \frac{\sin x}{x^{2}+1} $ בנקודה נתונה בצורה תכנותית. זה יראה כך:

public static double Derivative(double x)
{
    DualNumber number = new DualNumber(x, 1);

    DualNumber function =
        DualNumber.Sin(number) / (number*number + DualNumber.ONE);

    return function.DerivativeValue;
}

או אם נרצה לחשב את הנגזרת של הפונקציה $ e^{-\frac{x^{2}}{2}}\sin x $ זה יראה ככה:

public static double Derivative(double x)
{
    DualNumber number = new DualNumber(x, 1);

    DualNumber function =
        DualNumber.Exp(number * number * new DualNumber(-0.5, 0)) *
        DualNumber.Sin(number);

    return function.DerivativeValue;
}

כמובן, אפשר לשפר את הסינטקס, אבל זה ממש מגניב!

---

היתרון הגדול בשיטה היא הפשטות שלה: מאוד קל לגזור ככה פונקציות – פשוט נכתוב ביטוי מתמטי ונקבל חינם את החישוב של הנגזרת שלו.

יש מספר חסרונות לשיטה:

• אם מעניין אותנו נגזרת מסדר גבוה יותר, אפשר לעשות טריק דומה וליצור שלשות שבקואורדינטה השלישית שלהן מכבדות את חוקי הנגזרת השנייה. הבעיה שהחוקים עבור נגזרת שנייה ושלישית מתחילים להיות מאוד מסובכים ולא פשוטים כפי שהיה כאן.
• בנוסף, השיטה נותנת את הערך בנקודה, אבל לא נותנת לנו ביטוי מתמטי שמייצג את הנגזרת. החסרון הוא שלא נוכל למשל לפשט את הביטוי או לבצע בו אופטימיזציות כך שיהיה יותר מדויק/יותר מהיר וכו’.
  • בנוסף, לא נוכל להסתכל בעין מזוינת ולבצע בו ניתוח מתמטי – איפה הוא מוגדר, מתי הוא מקבל מינימום או מקסימום וכו’.
• לבסוף, אנחנו תלויים בדיוק של המחשב בו אנחנו עובדים – במקרים מסוימים אנחנו עשויים לסטות בהרבה מהתוצאה המדויקת של הנגזרת של הפונקציה בנקודה, וזאת רק מאחר ובצענו למשל יותר מדי פעולות חילוק או פעולות שלא שומרות על הדיוק.

שיהיה סופ"ש דואלי טוב.