CODEWITHSUNDEEP
pythonnumpymatrixphysicsnumerical-methods
Nam Kang
Nam Kang

Reputation: 3

I need help setting up matrices to solve using Gaussian elimination in Python

Apologies in advance to those who has to read through my poor coding skill

The objective of this coding is to first develop a 17x17 matrix and solve for the 17 unknowns using methods presented in linear algebra.

The part I am having the most difficulty is:

  1. implementing 2 counters i and j, where the value of i will increase once the value of j reaches its limit and goes back to 0 again.

  2. Lastly, being able to insert new values to a single array for later manipulation. I tried using np.insert, np.hstack, np.vstack, np.append, etc could not work it.

So i can generate matrix that looks like 

x11 x12 x13....x1j 
x21 .......... x2j
xi1............xij

here is some attempt 

import numpy as np
import math as mt
r=[2,2.8,3.2,3.5,3.7,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.7,3.5,3.2,2.8,2]
n=np.linspace(1,17,17)
m=np.linspace(1,17,17)
i=0
k=np.array([])
l=1
k2=[]
while i <=18:
    for j in range(17):
        h1=mt.sqrt(r[i]**2+(l*(n[i]-m[j])+l/2)**2)
        h2=mt.sqrt(r[i]**2+(l*(n[i]-m[j])-l/2)**2)
        h=h1-h2    
        k2.append(h)
        i=i+1

I am trying to obtain stokes' stream function in axially symmetrical flow for those who are interested,

I will appreciate any type of feedback, please guide me in the right direction

Upvotes: 0

Views: 583

Answers (2)

Cadmium
Cadmium

Reputation: 66

Your code suffers from two mistakes. The first that in Python, you start counting from zero; you may think of your matrix as having 17 rows, 1 to 17, but Python sees it as going from 0 to 16. The second is that when working with numpy, you should build your array first, and then insert your calculated values. There's a good explanation of why here:(How do I create an empty array/matrix in NumPy?).

I made r an array for consistency's sake, and I inserted the calculated values into k2. I'm not sure k was for. 

import numpy as np
import math as mt

r=np.array([2,2.8,3.2,3.5,3.7,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.7,3.5,3.2,2.8,2])
n=np.linspace(1,17,17)
m=np.linspace(1,17,17)
l=1

k2 = np.empty(shape=(17,17))
i=0
j=0
while i <=16:
    while j<=16: 
        h1=mt.sqrt(r[i]**2+(l*(n[i]-m[j])+l/2)**2)
        h2=mt.sqrt(r[i]**2+(l*(n[i]-m[j])-l/2)**2)
        h=np.array(h1-h2)    
        k2[i,j]= h
        j+=1
    j=0    
    i+=1

Upvotes: 1

Tonechas
Tonechas

Reputation: 13743

The code below is a vectorized solution to your problem:

import numpy as np

r = np.asarray([2,2.8,3.2,3.5,3.7,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.7,3.5,3.2,2.8,2])
l = 1

R = r.size
n, m = np.mgrid[1:R+1, 1:R+1]

h1 = np.sqrt(r[:, np.newaxis]**2 + (l*(n-m) + l/2.)**2)
h2 = np.sqrt(r[:, np.newaxis]**2 + (l*(n-m) - l/2.)**2)
k2 = h1 - h2

The result k2 is a 2-dimensional array rather than a vector:

>>> np.set_printoptions(precision=1)
>>> k2
array([[ 0. , -0.4, -0.7, -0.8, -0.9, -0.9, -0.9, -1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. ],
       [ 0.3,  0. , -0.3, -0.6, -0.7, -0.8, -0.9, -0.9, -0.9, -0.9, -1. , -1. , -1. , -1. , -1. , -1. , -1. ],
       [ 0.5,  0.3,  0. , -0.3, -0.5, -0.7, -0.8, -0.8, -0.9, -0.9, -0.9, -0.9, -1. , -1. , -1. , -1. , -1. ], 
       [ 0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.8, -0.8, -0.9, -0.9, -0.9, -0.9, -0.9, -1. , -1. , -1. ],
       [ 0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7, -0.8, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -1. ],
       [ 0.8,  0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7, -0.8, -0.8, -0.9, -0.9, -0.9, -0.9, -0.9],
       [ 0.8,  0.8,  0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7, -0.8, -0.8, -0.9, -0.9, -0.9, -0.9],
       [ 0.9,  0.8,  0.8,  0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7, -0.8, -0.8, -0.9, -0.9, -0.9],
       [ 0.9,  0.9,  0.8,  0.8,  0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7, -0.8, -0.8, -0.9, -0.9],
       [ 0.9,  0.9,  0.9,  0.8,  0.8,  0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7, -0.8, -0.8, -0.9],
       [ 0.9,  0.9,  0.9,  0.9,  0.8,  0.8,  0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7, -0.8, -0.8],
       [ 0.9,  0.9,  0.9,  0.9,  0.9,  0.8,  0.8,  0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7, -0.8],
       [ 1. ,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.8,  0.7,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6, -0.7],
       [ 1. ,  1. ,  1. ,  0.9,  0.9,  0.9,  0.9,  0.9,  0.8,  0.8,  0.6,  0.5,  0.3,  0. , -0.3, -0.5, -0.6],
       [ 1. ,  1. ,  1. ,  1. ,  1. ,  0.9,  0.9,  0.9,  0.9,  0.8,  0.8,  0.7,  0.5,  0.3,  0. , -0.3, -0.5],
       [ 1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  0.9,  0.9,  0.9,  0.9,  0.8,  0.7,  0.6,  0.3,  0. , -0.3],
       [ 1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  1. ,  0.9,  0.9,  0.9,  0.8,  0.7,  0.4,  0. ]])

Hopefully this is the result you were looking for.

Notice that in order to save space, only one decimal digit is displayed.

You may find it helpful to have a look on the description of the function mgrid and the object newaxis in Numpy's documentation to figure out how this code works.

Upvotes: 0

Related Questions