Sharedfirst-steps / MODEL.ipynbOpen in CoCalc
from sympy.abc import x, y, f, j
import sympy as sm

sm.init_printing()
sm.var('c1 c2 u0 u1 u2 u3 p0 p1 p2 p3 p11 p22 p33')  # This is the easiest way to define c as a standard symbol
\left ( c_{1}, \quad c_{2}, \quad u_{0}, \quad u_{1}, \quad u_{2}, \quad u_{3}, \quad p_{0}, \quad p_{1}, \quad p_{2}, \quad p_{3}, \quad p_{11}, \quad p_{22}, \quad p_{33}\right )
p0 = 1
p1 = 2
p2 = 0.5
p3 = -1
p11 = 0.1
p22 = 0.2
p33 = -0.05
x = p0*u0 + p1*u1 +p2*u2 +p3*u3 + p1*u1**2 +p2*u2**2 +p3*u3**2
x
u_{0} + 2 u_{1}^{2} + 2 u_{1} + 0.5 u_{2}^{2} + 0.5 u_{2} - u_{3}^{2} - u_{3}
c1 = 1.2
c2 = -0.8
y = c1*x + c2*f
y
- 0.8 f + 1.2 u_{0} + 2.4 u_{1}^{2} + 2.4 u_{1} + 0.6 u_{2}^{2} + 0.6 u_{2} - 1.2 u_{3}^{2} - 1.2 u_{3}
members = [u0, u1, u2, u3]
u = sm.FiniteSet(*members)
u_iterable = iter(u)
u = sm.IndexedBase('u')
p = sm.IndexedBase('p')
sum1 = sm.Sum(p[j]*u[j], (j, 0, 3))
sum2 = sm.Sum(p[j]*u[j]**2, (j, 1, 3))
x = sum1 + sum2
x
\sum_{j=0}^{3} p_{j} u_{j} + \sum_{j=1}^{3} p_{j} u_{j}^{2}
#p[0]=1
x.doit().evalf()
p_{0} u_{0} + p_{1} u_{1}^{2} + p_{1} u_{1} + p_{2} u_{2}^{2} + p_{2} u_{2} + p_{3} u_{3}^{2} + p_{3} u_{3}
u0 = 1
interv = []
interv.append(   sm.Interval(-5, 10)   )
interv.append(   sm.Interval(-7, 2)   )
interv.append(   sm.Interval(2, 13)   )
interv
\left [ \left[-5, 10\right], \quad \left[-7, 2\right], \quad \left[2, 13\right]\right ]
lboundaries = []
rboundaries = []
for i in range(3):
    u1min, u1max = interv[i].boundary
    lboundaries.append(u1min)
    rboundaries.append(u1max)
cons1 = interv[0].contains(u1); cons1
u_{1} \geq -5 \wedge u_{1} \leq 10
cons2 = interv[1].contains(u2); cons2
u_{2} \geq -7 \wedge u_{2} \leq 2
cons3 = interv[2].contains(u3); cons3
u_{3} \geq 2 \wedge u_{3} \leq 13
f = sm.Function('f')
sm.Sum(f(j), (j, 0, 3)).doit()
f{\left (0 \right )} + f{\left (1 \right )} + f{\left (2 \right )} + f{\left (3 \right )}
M = sm.Matrix([[1,10,2,13], [1,-5,2,2], [1,10,-7,2], [1,-5,-7,13]]); M
\left[\begin{matrix}1 & 10 & 2 & 13\\1 & -5 & 2 & 2\\1 & 10 & -7 & 2\\1 & -5 & -7 & 13\end{matrix}\right]
V = []
for i in range(4):
    for j in range(4):
        print(i,j)#if s: V[i,j] = lboundaries[i] else: rboundaries[i]
0 0 0 1 0 2 0 3 1 0 1 1 1 2 1 3 2 0 2 1 2 2 2 3 3 0 3 1 3 2 3 3
[  [lboundaries[i], rboundaries[i]] for i in range(3)  ]
\left [ \left [ -5, \quad 10\right ], \quad \left [ -7, \quad 2\right ], \quad \left [ 2, \quad 13\right ]\right ]
us = [lboundaries, rboundaries]; 
aa=[[row[i] for row in us] for i in range(3)]
aa
\left [ \left [ -5, \quad 10\right ], \quad \left [ -7, \quad 2\right ], \quad \left [ 2, \quad 13\right ]\right ]
VN = [[row[1-(n+i)%2] for i in range(4)] for n, row in enumerate(aa) ]
VN.insert(0, [1 for i in range(4)])
VN=[[row[i] for row in VN] for i in range(4)]
V1 = sm.Matrix(VN); V1
\left[\begin{matrix}1 & 10 & -7 & 13\\1 & -5 & 2 & 2\\1 & 10 & -7 & 13\\1 & -5 & 2 & 2\end{matrix}\right]
V = [[x,y] for x in lboundaries for y in rboundaries if x == y]
V1 = sm.Matrix(V); V1
[num for elem in us for num in elem]
u10 = (u1min + u1max)/2
u10.expand(basic=True)

sum( interv1.boundary ) / 2
from sympy.stats import P, E, variance, Normal, density, sample
f = Normal('f', 0, 1) 
sample(f)
variance(f)
density(f)(j)
from sympy.plotting import plot
pl1 = plot(density(f)(j), (j, -5, 5))
pl1

Exam

y1=[4.5, 3, 2, 0.5]
y2=[5.5, 5.5, 2, 1.5]
x0=[1, 1, 1, 1]
x1=[-1, 1, -1, 1]
x2=[-1, -1, 1, 1]
xs = [x0, x1, x2]
ysr=[]
for i in range(0,4):
    ysr.append( (y1[i]+y2[i])/2)
    
ysr
bj=[]
bj = [sum(    xs[j][i]*ysr[i] for i in range(4)   )/4 for j in range(3)]  
bj