Assignment 4 - Final Submission.ipynb

Kernel: Python 2 (SageMath)

שמות: רועי דביר

import sympy as sp
import math
from sympy.plotting import plot
sp.init_printing()
%matplotlib -- inline

עבודת סיכום SymPy

תרגיל ראשון

"יוחנן בוילר" הוא מהנדס בחברה המייצרת דודי מים חמים. במסגרת עבודתו נידרש יוחנן לתכנן דוד גלילי המיועד לחמום 150 ליטר (1 ליטר = 1000 סמ"ק) מים. יוחנן נידרש לתכנן את מימדי הדוד באופן כזה, שאובדן החום ממנו לסביבה יהיה הקטן ביותר האפשרי. ככל ששטח פני הדוד יהיה קטן יותר כך אובדן החום ממנו לסביבה יקטן. בעבודה זו עליכם למצוא מהם המימדים האידיאלים לדוד (רדיוס וגובה).

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

r,h=sp.symbols("r h")
S=2*sp.pi*r*(r+h)
display(S)

$$2 \pi r \left(h + r\right)$$

ב. נתון גליל שניפחו V ורדיוסו r. מצאו (בעזרת Sympy) ביטוי לגובהו והגדירו ביטוי לשטח הפנים כתלות בנפח והרדיוס.

h,r,V = sp.symbols('h r V')
eq1=sp.Eq(V,sp.pi*h*r**2)
hVal=sp.solve(eq1,[h])[0]
Snew=S.subs([(h,hVal)])
display(Snew)

$$2 \pi r \left(\frac{V}{\pi r^{2}} + r\right)$$

ג. שרטטו גרף של שטח פני דוד שניפחו 150 ליטר כתלות ברדיוס (הקפידו על עבודה נכונה עם היחידות).

sp.plot(Snew.subs([(V,150000)]),(r,0,100*1000),xlabel="Radius[cm]",ylabel="Surface[$cm^2$]")
#Only the positive part is true, because there is no negative radius or negative surface

<sympy.plotting.plot.Plot at 0x7f73fff1a910>

ד. עבור כל אחד מהנפחים 160L 40L,80L,120L,ו- 200L מצאו את רדיוס הגליל עבורו שטח הפנים מינימאלי.

Sdiff=Snew.diff(r)
display(Sdiff)
eq2=sp.Eq(Sdiff,0)
rVal=sp.solve(eq2,[r])[0]
rVal1=rVal.subs([(V,40*1000)])
rVal2=rVal.subs([(V,80*1000)])
rVal3=rVal.subs([(V,120*1000)])
rVal4=rVal.subs([(V,160*1000)])
rVal5=rVal.subs([(V,200*1000)])
display(rVal,rVal1.evalf(),rVal2.evalf(),rVal3.evalf(),rVal4.evalf(),rVal5.evalf())

$$2 \pi r \left(- \frac{2 V}{\pi r^{3}} + 1\right) + 2 \pi \left(\frac{V}{\pi r^{2}} + r\right)$$

$$\frac{2^{\frac{2}{3}} \sqrt[3]{V}}{2 \sqrt[3]{\pi}}$$

$$18.5336108963043$$

$$23.3508864988147$$

$$26.7300923514395$$

$$29.4202734335627$$

$$31.6920288377454$$

ה. עבור כל אחד מהרדיוסים שמצאתם בסעיף הקודם מצאו את גובה הגליל המתאים.

eq3=sp.Eq(S.subs([(r,rVal)]),Snew.subs([(r,rVal)]))
hVal=sp.solve(eq3,[h])[0]
hVal1=hVal.subs([(V,40*1000)])
hVal2=hVal.subs([(V,80*1000)])
hVal3=hVal.subs([(V,120*1000)])
hVal4=hVal.subs([(V,160*1000)])
hVal5=hVal.subs([(V,200*1000)])
display(hVal,hVal1.evalf(),hVal2.evalf(),hVal3.evalf(),hVal4.evalf(),hVal5.evalf())

$$\frac{2^{\frac{2}{3}} \sqrt[3]{V}}{\sqrt[3]{\pi}}$$

$$37.0672217926085$$

$$46.7017729976294$$

$$53.460184702879$$

$$58.8405468671253$$

$$63.3840576754908$$

ו. שרטטו גרף של גובה הגליל בו מתקבל מינימום שטח פנים כתלות ברדיוס. מה שיפוע הגרף? מה המסקנה?

#כפי שניתן לראות h=2r
eq4=sp.Eq(r,rVal)
vValInr=sp.solve(eq4,[V])[0]
display(eq4,vValInr)
sp.plot(hVal.subs([(V,vValInr)]),(r,0,20))
#משום כך השיפוע הינו 2

$$r = \frac{2^{\frac{2}{3}} \sqrt[3]{V}}{2 \sqrt[3]{\pi}}$$

$$2 \pi r^{3}$$

<sympy.plotting.plot.Plot at 0x7f73ff9f68d0>

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

z=sp.symbols("z")
eq5=sp.Eq(hVal.subs([(r,rVal)]).simplify(),z*rVal)
display(eq5)
slope=sp.solve(eq5,[z])[0]
slope

$$\frac{2^{\frac{2}{3}} \sqrt[3]{V}}{\sqrt[3]{\pi}} = \frac{2^{\frac{2}{3}} \sqrt[3]{V} z}{2 \sqrt[3]{\pi}}$$

$$2$$

---

תרגיל שני

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

a,x=sp.symbols('a,x')
f=sp.exp(a*x)/x
df=f.diff(x)
f,df

$$\left ( \frac{e^{a x}}{x}, \quad \frac{a}{x} e^{a x} - \frac{e^{a x}}{x^{2}}\right )$$

g=sp.lambdify((x,a),df)
g(0.5,2),g(0.5,1)

$$\left ( 0.0, \quad -3.2974425414\right )$$

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

בהצלחה!

import sympy as sp
x=sp.symbols('x')
f = sp.cos(x)*sp.sin(x)+2*sp.sin(x)
sp.plot(f,(x,-10,10),xlabel="x",ylabel="y",legend=True,axis=True,aspect_ratio=('auto'),autoscale=True,axis_center=('auto'))
def GetDiff(f):
    return sp.diff(f,x)
def func(f):
    PointKitzon = []
    PointMin = []
    PointMax = []
    PointPitulMin = []
    PointPitulMax = []
    yHatzava = 0
    count = 0
    fDf = GetDiff(f)
    sDf = GetDiff(fDf)
    tDf = GetDiff(sDf)
    eq1 = sp.Eq(fDf,0)
    sol1=sp.solve(eq1,[x])
    i=0
    for item in sol1:
        if(item.is_complex):
            del sol1[i]
            i-=1
        i+=1
    eq2 = sp.Eq(sDf,0)
    sol2=sp.solve(eq2,[x])
    for item in sol2:
        if(item.is_complex):
            del sol2[i]
            i-=1
        i+=1
    display(sol1,sol2)
    if(sol1):
        for i in range(len(sol1)):
            x1=sol1[i]
            yHatzava = sDf.subs({x:x1})
            if(yHatzava > 0):
                PointMin.append([x1,f.subs({x:x1})])
                PointKitzon.append([x1,f.subs({x:x1})])
            elif(yHatzava < 0):
                PointMax.append([x1,f.subs({x:x1})])
                PointKitzon.append([x1,f.subs({x:x1})])
            for i in range(len(PointMax)):
                x1=PointMax[i]
                print("Max: ",x1[0].simplify().evalf(),x1[1].simplify().evalf())
            for i in range(len(PointMin)):
                x1=PointMin[i]
                print("Min: ",x1[0].simplify().evalf(),x1[1].simplify().evalf())
    else:
        print("There are no max and min points.")
        count+=1
    if(sol2):
        for i in range(len(sol2)):
            x1=sol2[i]
            yHatzava = tDf.subs({x:x1})
            if(yHatzava > 0):
                PointPitulMax.append([x1,f.subs({x:x1})])
                PointKitzon.append([x1,f.subs({x:x1})])
            elif(yHatzava < 0):
                PointPitulMin.append([x1,f.subs({x:x1})])
                PointKitzon.append([x1,f.subs({x:x1})])
        for i in range(len(PointPitulMax)):
            x1=PointPitulMax[i]
            print("Pitul Max: ",x1[0].simplify().evalf(),x1[1].simplify().evalf())
        for i in range(len(PointPitulMin)):
            x1=PointPitulMin[i]
            print("Pitul Min: ",x1[0].simplify().evalf(),x1[1].simplify().evalf())
    else:
        print("There are no pitul points.")
        count+=1
    if(count == 2):
        print("There are no kitzon points.")
    return(PointKitzon)
try:
    func(f)
except TypeError:
    print("The function is invalid or there was a problem")

[-2*atan(sqrt(-3 + 2*sqrt(3))), 2*atan(sqrt(-3 + 2*sqrt(3)))]

[2*pi/3]

('Min: ', -1.19606189408616, -2.20183473752081)
('Max: ', 1.19606189408616, 2.20183473752081)
('Min: ', -1.19606189408616, -2.20183473752081)
('Pitul Max: ', 2.09439510239320, 1.29903810567666)