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