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
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
c1 = 1.2
c2 = -0.8
y = c1*x + c2*f
y
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
#p[0]=1
x.doit().evalf()
u0 = 1
interv = []
interv.append(   sm.Interval(-5, 10)   )
interv.append(   sm.Interval(-7, 2)   )
interv.append(   sm.Interval(2, 13)   )
interv
lboundaries = []
rboundaries = []
for i in range(3):
u1min, u1max = interv[i].boundary
lboundaries.append(u1min)
rboundaries.append(u1max)
cons1 = interv[0].contains(u1); cons1
cons2 = interv[1].contains(u2); cons2
cons3 = interv[2].contains(u3); cons3
f = sm.Function('f')
sm.Sum(f(j), (j, 0, 3)).doit()
M = sm.Matrix([[1,10,2,13], [1,-5,2,2], [1,10,-7,2], [1,-5,-7,13]]); M
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)  ]
us = [lboundaries, rboundaries];
aa=[[row[i] for row in us] for i in range(3)]
aa
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
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