Difference between revisions of "Metode Numerik - Josiah Enrico S (1906356286)"

From ccitonlinewiki
Jump to: navigation, search
(Tugas Besar Metode Numerik - Metode Numerik/23 Desember 2020)
(Jawaban UAS 2021 - Josiah Enrico S)
 
(66 intermediate revisions by the same user not shown)
Line 1,000: Line 1,000:
  
  
== Tugas Besar Metode Numerik - Metode Numerik/23 Desember 2020 ==
+
== [[Tugas Besar Metode Numerik (Latihan)- Metode Numerik/23 Desember 2020]] ==
  
 +
==[[Tugas Besar Metode Numerik - Josiah Enrico S]]==
  
'''Objektif:'''
 
  
- Mengoptimasi harga pembuatan rangka truss sederhana dengan memvariasi dimensi dan elastisitas material.
+
==Jawaban UAS 2021 - Josiah Enrico S==
 +
<gallery mode="slideshow">
 +
UAS Metnum02 Josiah Enrico S-01.jpg
 +
UAS Metnum02 Josiah Enrico S-02.jpg
 +
UAS Metnum02 Josiah Enrico S-03.jpg
 +
UAS Metnum02 Josiah Enrico S-04.jpg
 +
UAS Metnum02 Josiah Enrico S-05.jpg
 +
UAS Metnum02 Josiah Enrico S-06.jpg
 +
UAS Metnum02 Josiah Enrico S-07.jpg
 +
UAS Metnum02 Josiah Enrico S-08.jpg
 +
UAS Metnum02 Josiah Enrico S-09.jpg
 +
UAS Metnum02 Josiah Enrico S-10.jpg
 +
</gallery>
  
'''Geometri dan Load'''
+
Revisi dan Narasi
 +
=== No 1 ===
  
[[File:Tugas Besar Metnum Geometri Jos.jpg|center]]
+
Langkah-langkah optimasi
  
'''Constraint:'''
+
# Tentukan tujuan optimasi lalu definisikan semua parameter yang ada.
 +
# Menentukan boundary dan constraint pada struktur
 +
# Menentukan asumsi dalam permodelan struktur
 +
# Melakukan analisa struktur (komputasi)
 +
# Setelah mendapat persamaan utama, variasikan nilai objektif yang akan dioptimasi
 +
# Gunakan metode optimasi yang sesuai
 +
# Setelah mendapatakan nilai optimasi, lakukan pengujian variabel dengan analisa
 +
# Jika hasilnya memuaskan, optimasi selesai. Jika hasilnya belum memuaskan, ulangi dari langkah 6
  
- Spesifikasi L (Panjang) dan geometri rangka truss
+
=== No 2 ===
  
- Gaya beban terhadap struktur (1000 N dan 2000 N)
+
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
 +
|-
 +
| '''Tujuan Optimasi:'''
 +
|
 +
* Mengoptimasi kekuatan material struktur dan dinding water tank
 +
|-
 +
| '''Asumsi:'''
 +
|
 +
* Fluida Statis
 +
* Gravitasi: 9,81 m/s2
 +
* Fluida: air murni (rho = 1000 kg/m3)
 +
|-
 +
| '''Dalil Fisika:'''
 +
|
 +
* Gaya Hidrostatis F = rho*g*V
 +
* Stress = E*Strain
 +
|}
  
'''Asumsi:'''
+
=== No 3 ===
  
- Variasi Stiffness terikat dengan variabel area. Memvariasikan Elastisitas tergolong sulit karena setiap material memiliki range yang tidak teratur dan dalam satu material yang sejenis (struktur biaya tetap) tidak terjadi perubahan nilai elastisitas yang berbanding lurus dengan perubahan biaya.
+
[[File:UAS Metnum Water Tower.jpg|thumb|600px|right|'''Nomor 3''' Ilustrasi struktur truss water tower (3 elemen)]]
 +
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
 +
|+ Matriks Persamaan Aljabar
 +
|-
 +
| Persamaan aljabar lokal tiap elemen
 +
| Persamaan aljabar global tiap elemen
 +
|-
 +
| style='border-style: none right none left;' |
  
- Beban akan terdistribusi hanya pada point penghubung (karena bersifat truss)
+
Truss A (1 to 2)
  
- Safety factor bernilai 2.
+
            [ 0  0  0  0 ] (1x)
 +
k = (A*E/L) [ 0  1  0  1 ] (1y)
 +
            [ 0  0  0  0 ] (2x)
 +
            [ 0 -1  0 -1 ] (2y)
  
- Batas displacement 0,001 m sebelum buckling(pada kolom paling atas)
+
| style='border-style: none right none left;' |
  
Koleksi Data:
+
Truss A (1 to 2)
  
 +
            [ 0  0  0  0  0  0  0  0 ] (1x)
 +
            [ 0  1  0  1  0  0  0  0 ] (1y)
 +
            [ 0  0  0  0  0  0  0  0 ] (2x)
 +
k = (A*E/L) [ 0 -1  0 -1  0  0  0  0 ] (2y)
 +
            [ 0  0  0  0  0  0  0  0 ] (3x)
 +
            [ 0  0  0  0  0  0  0  0 ] (3y)
 +
            [ 0  0  0  0  0  0  0  0 ] (4x)
 +
            [ 0  0  0  0  0  0  0  0 ] (4y)
 +
|-
 +
| style='border-style: none right none left;' |
  
[[File:Tugas Besar Metnum Dimension Jos.jpg]] [[File:Tugas Besar Metnum Model Jos.jpg]]
+
Truss B (2 to 3)
  
[[File:Tugas Besar Metnum AverageCostSS304 Jos.jpg]]
+
            [ 0.5  0.5 -0.5 -0.5 ] (2x)
 +
k = (A*E/L) [ 0.5  0.5 -0.5 -0.5 ] (2y)
 +
            [-0.5 -0.5  0.5  0.5 ] (3x)
 +
            [-0.5 -0.5  0.5  0.5 ] (3y)
  
[[File:Tugas Besar Metnum AverageCostASTMA36 Jos.jpg]]
+
| style='border-style: none right none left;' |
  
 +
Truss B (2 to 3)
  
 +
            [ 0  0  0    0    0    0    0  0 ] (1x)
 +
            [ 0  0  0    0    0    0    0  0 ] (1y)
 +
            [ 0  0  0.5  0.5 -0.5 -0.5  0  0 ] (2x)
 +
k = (A*E/L) [ 0  0  0.5  0.5 -0.5 -0.5  0  0 ] (2y)
 +
            [ 0  0 -0.5 -0.5  0.5  0.5  0  0 ] (3x)
 +
            [ 0  0 -0.5 -0.5  0.5  0.5  0  0 ] (3y)
 +
            [ 0  0  0    0    0    0    0  0 ] (4x)
 +
            [ 0  0  0    0    0    0    0  0 ] (4y)
 +
|-
 +
| style='border-style: none right none left;' |
  
Resources:
+
Truss C (2 to 4)
  
http://www.matweb.com/search/DataSheet.aspx?MatGUID=1748ca73d11e4353b2aa700bfb119dfb
+
            [ 0.5 -0.5 -0.5  0.5 ] (2x)
 +
k = (A*E/L) [-0.5  0.5  0.5 -0.5 ] (2y)
 +
            [-0.5  0.5  0.5 -0.5 ] (4x)
 +
            [ 0.5 -0.5 -0.5  0.5 ] (4y)
  
http://www.matweb.com/search/datasheet.aspx?matguid=d1844977c5c8440cb9a3a967f8909c3a
+
| style='border-style: none right none left;' |
  
https://mitarcahyaabadai.wordpress.com/daftar-harga-besi-siku-2018/
+
Truss C (2 to 4)
  
https://duniabahanbangunanbandung.blogspot.com/p/harga-besi-siku-stainless-steel.html
+
            [ 0  0  0    0    0  0  0    0  ] (1x)
 +
            [ 0  0  0    0    0  0  0    0  ] (1y)
 +
            [ 0  0  0.5 -0.5  0  0 -0.5  0.5 ] (2x)
 +
k = (A*E/L) [ 0  0 -0.5  0.5  0  0  0.5 -0.5 ] (2y)
 +
            [ 0  0  0    0    0  0  0    0  ] (3x)
 +
            [ 0  0  0    0    0  0  0    0  ] (3y)
 +
            [ 0  0 -0.5  0.5  0  0  0.5 -0.5 ] (4x)
 +
            [ 0  0  0.5 -0.5  0  0 -0.5  0.5 ] (4y)
  
http://www.wermac.org/steel/dim_angle_eq.html
+
|}
 +
 
 +
=== No 4 ===
 +
 
 +
<gallery mode=packed heights=500px>
 +
Kuis Metnum1 Jos.jpeg|'''Nomor 4''' Flowchart Analisa Struktur
 +
Verifikasi Metnum Jos.jpeg|'''Nomor 4''' Verifikasi Gaya Reaksi
 +
</gallery>
 +
 
 +
=== No 5 ===
 +
 
 +
Fungsi Objektif
 +
 
 +
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
 +
|-
 +
| displacement=
 +
| [A*E/L][Uglobal]* '''deltaL''' = [F]
 +
|-
 +
| Reaction Force=
 +
| '''R''' = [A*E/L][Uglobal]*deltaL - [F]
 +
|-
 +
| Stress=
 +
| '''Stress''' = E*(deltaL)/A
 +
|}
 +
 
 +
Constraint/Boundary
 +
 
 +
*Berdasarkan struktur di nomor 3, boundary condition terdapat pada titik 3 dan 4
 +
 
 +
=== No 6 ===
 +
 
 +
Asumsi:
 +
 
 +
*Area: 0.16 (Persegi 0.4 x 0.4 meter)
 +
*Elastisitas: 193 GPa (SS 304)
 +
*Yield: 215 MPa (SS 304)
 +
*rho air murni: 1000 kg/m3
 +
*Gravitasi: 9,81 m/s2
 +
*Volume air: 200000 gallon / 1136 m3
 +
 
 +
=== No 7 ===
 +
 
 +
<gallery mode=packed heights=300px>
 +
UAS Metnum Result.jpg|'''Nomor 7''' Hasil Simulasi Modelica
 +
</gallery>
 +
 
 +
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
 +
|-
 +
|
 +
'''Kode Modelica'''
 +
 
 +
model UAS_METNUM
 +
 +
//define initial variable
 +
parameter Integer Points=4; //Number of Points
 +
parameter Integer Trusses=3; //Number of Trusses
 +
parameter Real Yield=215e6; //Yield Strength (Pa)
 +
parameter Real Area=0.16;  //Area Block Profile (0.4 * 0.4) (m2)
 +
parameter Real Elas=193e9;    //Elasticity SS 304  (Pa)
 +
parameter Real Volume=1135; //m3 from 200000 gallons
 +
parameter Real rho=1000; //kg/m3
 +
 +
parameter Real W=Volume*rho*9.81;
 +
 +
//define connection
 +
parameter Integer C[Trusses,2]=[1,2;
 +
                                2,3;
 +
                                2,4];
 +
                                                             
 +
//define coordinates (please put orderly)
 +
parameter Real P[Points,3]=[0,37,0;
 +
                            0,8,0;
 +
                            -8,0,0;
 +
                            8,0,0];
 +
                           
 +
//define external force (please put orderly)
 +
parameter Real F[Points*3]={0,-W,0,
 +
                            0,0,0,
 +
                            0,0,0,
 +
                            0,0,0};
 +
 +
//define boundary
 +
parameter Integer b[:]={3,4};
 +
 +
//solution
 +
Real displacement[N], reaction[N];
 +
Real check[3];
 +
 +
Real stress1[Trusses];
 +
Real safety[Trusses];
 +
Real dis[3];
 +
Real Str[3];
 +
 +
protected
 +
parameter Integer N=3*Points;
 +
Integer boundary[3*size(b,1)]=cat(1,(3*b).-2,(3*b).-1,3*b);
 +
Real q1[3], q2[3], g[N,N], G[N,N], G_star[N,N], id[N,N]=identity(N), cx, cy, cz, L, X[3,3];
 +
Real err=10e-10, ers=10e-4;
 +
 +
algorithm
 +
//Creating Global Matrix
 +
G:=id;
 +
for i in 1:Trusses loop
 +
for j in 1:3 loop
 +
  q1[j]:=P[C[i,1],j];
 +
  q2[j]:=P[C[i,2],j];
 +
end for;
 +
     
 +
    //Solving Matrix
 +
    L:=Modelica.Math.Vectors.length(q2-q1);
 +
    cx:=(q2[1]-q1[1])/L;
 +
    cy:=(q2[2]-q1[2])/L;
 +
    cz:=(q2[3]-q1[3])/L;
 +
    X:=(Area*Elas/L)*[cx^2,cx*cy,cx*cz;
 +
                      cy*cx,cy^2,cy*cz;
 +
                      cz*cx,cz*cy,cz^2];
 +
 +
    //Transforming to global matrix
 +
    g:=zeros(N,N);
 +
    for m,n in 1:3 loop
 +
      g[3*(C[i,1]-1)+m,3*(C[i,1]-1)+n]:=X[m,n];
 +
      g[3*(C[i,2]-1)+m,3*(C[i,2]-1)+n]:=X[m,n];
 +
      g[3*(C[i,2]-1)+m,3*(C[i,1]-1)+n]:=-X[m,n];
 +
      g[3*(C[i,1]-1)+m,3*(C[i,2]-1)+n]:=-X[m,n];
 +
    end for; 
 +
 +
G_star:=G+g;
 +
G:=G_star;
 +
end for;
 +
 +
//Implementing boundary
 +
for i in boundary loop
 +
for j in 1:N loop
 +
  G[i,j]:=id[i,j];
 +
end for;
 +
end for;
 +
 +
//Solving displacement
 +
displacement:=Modelica.Math.Matrices.solve(G,F);
 +
 +
//Solving reaction
 +
reaction:=(G_star*displacement)-F;
 +
 +
//Eliminating float error
 +
for i in 1:N loop
 +
reaction[i]:=if abs(reaction[i])<=err then 0 else reaction[i];
 +
displacement[i]:=if abs(displacement[i])<=err then 0 else displacement[i];
 +
end for;
 +
 +
//Checking Force
 +
check[1]:=sum({reaction[i] for i in (1:3:(N-2))})+sum({F[i] for i in (1:3:(N-2))});
 +
check[2]:=sum({reaction[i] for i in (2:3:(N-1))})+sum({F[i] for i in (2:3:(N-1))});
 +
check[3]:=sum({reaction[i] for i in (3:3:N)})+sum({F[i] for i in (3:3:N)});
 +
 
 +
for i in 1:3 loop
 +
  check[i] := if abs(check[i])<=ers then 0 else check[i];
 +
end for;
 +
 +
//Calculating stress in each truss
 +
for i in 1:Trusses loop
 +
for j in 1:3 loop
 +
  q1[j]:=P[C[i,1],j];
 +
  q2[j]:=P[C[i,2],j];
 +
  dis[j]:=abs(displacement[3*(C[i,1]-1)+j]-displacement[3*(C[i,2]-1)+j]);
 +
end for;
 +
     
 +
    //Solving Matrix
 +
    L:=Modelica.Math.Vectors.length(q2-q1);
 +
    cx:=(q2[1]-q1[1])/L;
 +
    cy:=(q2[2]-q1[2])/L;
 +
    cz:=(q2[3]-q1[3])/L;
 +
    X:=(Elas/L)*[cx^2,cx*cy,cx*cz;
 +
                cy*cx,cy^2,cy*cz;
 +
                cz*cx,cz*cy,cz^2];
 +
   
 +
    Str:=(X*dis);
 +
    stress1[i]:=Modelica.Math.Vectors.length(Str);
 +
end for;
 +
 +
//Safety factor
 +
for i in 1:Trusses loop
 +
  if stress1[i]>0 then
 +
    safety[i]:=Yield/stress1[i];
 +
  else
 +
    safety[i]:=0;
 +
  end if;
 +
end for;
 +
 +
end UAS_METNUM;
 +
|}

Latest revision as of 21:24, 14 January 2021

Intro - Metode Numerik/11 November 2020

Tujuan belajar metode numerik yaitu agar penulis dapat:

1) Memahami konsep dan prinsip dasar metode numerik.

2) Memahami aplikasi metode numerik

3) Menyelesaikan persoalan teknik dengan metoda numerik

4) Mendapat nilai tambah sebagai manusia yang lebih beradab


Aplikasi Modelica - Metode Numerik/18 November 2020

Berikut ini adalah contoh penerapan aplikasi OpenModelica untuk membuat 4 algoritma metode numerik dalam mencari roots of equation (akar persamaan) dari:

f(x) = exp^(-x)-(x)

f'(x) = -exp^(-x)-1

error maksimum = 0.0000001

1) Newton Raphson (Terbuka)

model Newton_Raphson_Algorithm

parameter Real g=1; //guess
parameter Integer N=20; //max iteration
parameter Real er=0.0000001; //error maximum
Real a[N]; 
Real y[N];//function
Real ER[N]; //error
Real sol; //solution

algorithm

a[1]:=g;
y[1]:=a[1]-(exp(-a[1])-a[1])/(-exp(-a[1])-1);
ER[1]:=abs(1-a[1]/y[1]);

for i in 2:N loop
  a[i]:=y[i-1];
  y[i]:=a[i]-(exp(-a[i])-a[i])/(-exp(-a[i])-1);
  ER[i]:=abs(1-y[i-1]/y[i]);

  if ER[i]<er then
 sol:=y[i];
 break;
 end if;
end for;

end Newton_Raphson_Algorithm;

2) Secant (Terbuka)

model Secant_Algorithm

parameter Real a=0; //guess
parameter Real b=1; //guess
parameter Integer N=10; //max iteration
parameter Real er=0.0000001; //error maximum
Real A[N];
Real B[N];
Real y[N];
Real ER[N];
Real sol; //solution

algorithm

A[1]:=a;
B[1]:=b;
y[1]:=B[1]-(exp(-B[1])-B[1])*(A[1]-B[1])/((exp(-A[1])-A[1])-(exp(-B[1])-B[1]));
ER[1]:=abs(1-B[1]/y[1]);

for i in 2:N loop
 A[i]:=B[i-1];
 B[i]:=y[i-1];
 y[i]:=B[i]-(exp(-B[i])-B[i])*(A[i]-B[i])/((exp(-A[i])-A[i])-(exp(-B[i])-B[i]));
 ER[i]:=abs(1-y[i-1]/y[i]);

 if ER[i]<er then
 sol:=y[i];
 break;
 
 end if;
end for;

end Secant_Algorithm;

3) Bisection (Tertutup)

model Bisection_Algorithm

parameter Real a=0; //guess bawah
parameter Real b=1; //guess atas
parameter Integer N=50; //max iteration
parameter Real er=0.0000001; //error maximum
Real fa=(exp(-a)-a);
Real fb=(exp(-b)-b);
Real A[N];
Real B[N];
Real fy[N];
Real y[N];
Real ER[N];
Real sol; //solution

algorithm

if fa*fb<0 then

A[1]:=a;
B[1]:=b;
y[1]:=(A[1]+B[1])/2;
fy[1]:=exp(-y[1])-y[1];
ER[1]:=1;

for i in 2:N loop
 if fy[i-1]>0 then
 A[i]:=y[i-1];
 B[i]:=B[i-1];
 else
 A[i]:=A[i-1];
 B[i]:=y[i-1];
 end if;
   
 y[i]:=(A[i]+ B[i])/2;
 fy[i]:=exp(-y[i])-y[i];
 ER[i]:=abs(1-y[i-1]/y[i]);
 
 if ER[i]<er then
 sol:=y[i];
 break;
 end if;

end for;
end if;

end Bisection_Algorithm;

4) Regula Falsi (Tertutup)

model Regula_Falsi_Algorithm

parameter Real a=0; //guess bawah
parameter Real b=1; //guess atas
parameter Integer N=20; //max iteration
parameter Real er=0.0000001; //error maximum
Real A[N];
Real B[N];
Real fa[N];
Real fb[N];
Real fy[N];
Real y[N];
Real ER[N];
Real sol; //solution

algorithm

A[1]:=a;
B[1]:=b;
fa[1]:=exp(-A[1])-A[1];
fb[1]:=exp(-B[1])-B[1]; 

if fa[1]*fb[1]<0 then 

y[1]:=(A[1]*fb[1]-B[1]*fa[1])/(fb[1]-fa[1]);
fy[1]:=exp(-y[1])-y[1];
ER[1]:=1; 

for i in 2:N loop
 if fy[i-1]>0 then
 A[i]:=y[i-1];
 B[i]:=B[i-1];
 else
 A[i]:=A[i-1];
 B[i]:=y[i-1];
 end if;
 
 fa[i]:=exp(-A[i])-A[i];
 fb[i]:=exp(-B[i])-B[i];
 y[i]:=(A[i]*fb[i]-B[i]*fa[i])/(fb[i]-fa[i]);
 fy[i]:=exp(-y[i])-y[i];
 ER[i]:=abs(1-y[i-1]/y[i]);
 
 if ER[i]<er then
 sol:=y[i];
 break;
 end if;
end for;
end if;
 
end Regula_Falsi_Algorithm;

Berikut link Youtube yang berisi penjelasan tentang algoritma dan penulisan kode OpenModelica keempat metode diatas:

Fungsi Panggil dalam Modelica - Metode Numerik/18 November 2020

Fungsi panggil dalam aplikasi Modelica adalah metode membuat fungsi numerik dalam kelas function yang akan digunakan di dalam permodelan numerik utama (tipe 'class'). metode ini biasanya dipakai untuk menyederhanakan persoalan matematika yang kompleks, seperti persamaan aljabar simultan, sehingga lebih mudah diselesaikan. Dalam artikel ini akan ditampilkan 3 contoh persoalan numerik dengan solusi array atau lebih dari satu anggota:

1) Root of Equation (Newton Raphson dan Regula Falsi)

Grafik Modelica

Newton Raphson

...
algorithm

a[1]:=g;
y[1]:=a[1]-(Fungsi(a[1]))/(Der_Fungsi(a[1]));
ER[1]:=abs(1-a[1]/y[1]);

for i in 2:N loop
 a[i]:=y[i-1];
 y[i]:=a[i]-(Fungsi(a[i]))/(Der_Fungsi(a[1]));
 ER[i]:=abs(1-y[i-1]/y[i]);

 if ER[i]<er then
 sol:=y[i];
 break;
 end if;

Regula Falsi

...
fa[1]:=Fungsi(A[1]);
fb[1]:=Fungsi(B[1]);

if fa[1]*fb[1]<0 then

y[1]:=(A[1]*fb[1]-B[1]*fa[1])/(fb[1]-fa[1]);
fy[1]:=Fungsi(y[1]);
ER[1]:=1;
...
 
 fa[i]:=Fungsi(A[i]);
 fb[i]:=Fungsi(B[i]);
 y[i]:=(A[i]*fb[i]-B[i]*fa[i])/(fb[i]-fa[i]);
 fy[i]:=Fungsi(y[i]);
 ER[i]:=abs(1-y[i-1]/y[i]);
...

fungsi

function Fungsi
input Real X;
output Real Y;

algorithm
Y:=exp(-X)-X;
end Fungsi;

fungsi turunan

function Der_Fungsi
input Real X;
output Real Y;

algorithm
Y:=-exp(-X)-1;
end Der_Fungsi;

2) Heat Diffusion

Referensi: Versteeg, H. K., Malalasekera, W. (2007). An Introduction to Computational Fluid Dynamics. 2nd Edition. Harlow: Pearson (Example 4.2)

Grafik Modelica

Persamaan

class TesAljabar
parameter Real A[5,5]=[375,-125,0,0,0;
                    -125,250,-125,0,0;
                    0,-125,250,-125,0;
                    0,0,-125,250,-125;
                    0,0,0,-125,375];

parameter Real B[5]={29000,4000,4000,4000,54000};
Real X[5];

equation
X=matriks5(A,B);
end TesAljabar;

fungsi

function matriks5
input Real A[5,5];
input Real B[5];
output Real X[5];

algorithm
X:=Modelica.Math.Matrices.solve(A,B);
end matriks5;

3) Combined Spring Stiffness

Referensi: http://www.sharetechnote.com/html/EngMath_Matrix_FEM_ComplicatedSpring.html

Grafik Modelica

Persamaan

class Latihan_Matriks

parameter Real c1=1; //spring 1
parameter Real c2=2; //spring 2
parameter Real c3=3; //spring 3
parameter Real c4=4; //spring 4
parameter Real c5=5; //spring 5
parameter Real k1[3,3]=c1*[1,0,0;
                           0,0,0;
                           0,0,0];
parameter Real k2[3,3]=c2*[1,-1,0;
                           -1,1,0;
                           0,0,0];
parameter Real k3[3,3]=c3*[1,-1,0;
                           -1,1,0;
                           0,0,0];
parameter Real k4[3,3]=c4*[1,0,-1;
                           0,0,1;
                           -1,0,1];
parameter Real k5[3,3]=c5*[0,0,0;
                           0,1,-1;
                           0,-1,1];
  
parameter Real B[3]={100,200,300};
Real K[3,3]; 
Real s[3];
 
algorithm
K:=k1+k2+k3+k4+k5;
s:=matriks3(K,B);

end Latihan_Matriks;

fungsi

function matriks3
input Real A[3,3];
input Real B[3];
output Real X[3];

algorithm
X:=Modelica.Math.Matrices.solve(A,B);
end matriks3;

Berikut link Youtube yang berisi penjelasan tentang algoritma dan penulisan kode OpenModelica latihan diatas:

Berikut juga adalah pemodelan Modelica untuk metode Naive Gauss:

Finite Element Method for Trusses - Metode Numerik/2 Desember 2020

Berikut 2 contoh penggunaan aplikasi Modelica untuk menyelesaikan perhitungan displacement dan reaction force pada trusses:

Trusses Problem 1 (Example 3.1)

Soal Trusses 1 Jos.jpg
Grafik Displacement
Grafik Reaction Forces

Persamaan

model Trusses

parameter Integer N=10; //Global matrice = 2*points connected
parameter Real A=8;
parameter Real E=1.9e6;
Real G[N,N]; //global
Real Ginitial[N,N]; //global
Real Sol[N]; //global dispplacement
Real X[N]={0,0,0,0,0,0,0,-500,0,-500};
Real R[N]; //global reaction force
Real SolMat[N,1];
Real XMat[N,1];

//boundary coundition
Integer b1=1;
Integer b2=3;

//truss 1
parameter Real X1=0; //degree between truss
Real k1=A*E/36;
Real K1[4,4]; //stiffness matrice
Integer p1a=1;
Integer p1b=2;
Real G1[N,N];

//truss 2
parameter Real X2=135; //degree between truss
Real k2=A*E/50.912;
Real K2[4,4]; //stiffness matrice
Integer p2a=2;
Integer p2b=3;
Real G2[N,N];

//truss 3
parameter Real X3=0; //degree between truss
Real k3=A*E/36;
Real K3[4,4]; //stiffness matrice
Integer p3a=3;
Integer p3b=4;
Real G3[N,N];

//truss 4
parameter Real X4=90; //degree between truss
Real k4=A*E/36;
Real K4[4,4]; //stiffness matrice
Integer p4a=2;
Integer p4b=4;
Real G4[N,N];

//truss 5
parameter Real X5=45; //degree between truss
Real k5=A*E/50.912;
Real K5[4,4]; //stiffness matrice
Integer p5a=2;
Integer p5b=5;
Real G5[N,N];

//truss 6
parameter Real X6=0; //degree between truss
Real k6=A*E/36;
Real K6[4,4]; //stiffness matrice
Integer p6a=4;
Integer p6b=5;
Real G6[N,N];

/*
for each truss, please ensure pXa is lower then pXb (X represents truss element number)
*/

algorithm

//creating global matrice
K1:=Stiffness_Matrices(X1);
G1:=k1*Local_Global(K1,N,p1a,p1b);

K2:=Stiffness_Matrices(X2);
G2:=k2*Local_Global(K2,N,p2a,p2b);

K3:=Stiffness_Matrices(X3);
G3:=k3*Local_Global(K3,N,p3a,p3b);

K4:=Stiffness_Matrices(X4);
G4:=k4*Local_Global(K4,N,p4a,p4b);

K5:=Stiffness_Matrices(X5);
G5:=k5*Local_Global(K5,N,p5a,p5b);

K6:=Stiffness_Matrices(X6);
G6:=k6*Local_Global(K6,N,p6a,p6b);

G:=G1+G2+G3+G4+G5+G6;
Ginitial:=G;

//implementing boundary condition
for i in 1:N loop
 G[2*b1-1,i]:=0;
 G[2*b1,i]:=0;
 G[2*b2-1,i]:=0;
 G[2*b2,i]:=0;
end for;

G[2*b1-1,2*b1-1]:=1;
G[2*b1,2*b1]:=1;
G[2*b2-1,2*b2-1]:=1;
G[2*b2,2*b2]:=1;

//solving displacement
Sol:=Gauss_Jordan(N,G,X);

//solving reaction force
SolMat:=matrix(Sol);
XMat:=matrix(X);
R:=Reaction_Trusses(N,Ginitial,SolMat,XMat);

end Trusses;


Trusses Problem 2 (Homework)

Soal Trusses 2 Jos.jpg
Grafik Displacement
Grafik Reaction Forces

Persamaan

class Trusses_HW

parameter Integer N=8; //Global matrice = 2*points connected
parameter Real A=0.001; //Area m2
parameter Real E=200e9; //Pa
Real G[N,N]; //global
Real Ginitial[N,N]; //global
Real Sol[N]; //global dispplacement
Real X[N]={0,0,-1035.2762,-3863.7033,0,0,-1035.2762,-3863.7033};
Real R[N]; //global reaction force
Real SolMat[N,1];
Real XMat[N,1];

//boundary condition
Integer b1=1;
Integer b2=3;

//truss 1
parameter Real X1=0; //degree between truss
Real k1=A*E/1;
Real K1[4,4]; //stiffness matrice
Integer p1a=1;
Integer p1b=2;
Real G1[N,N];

//truss 2
parameter Real X2=0; //degree between truss
Real k2=A*E/1;
Real K2[4,4]; //stiffness matrice
Integer p2a=2;
Integer p2b=3;
Real G2[N,N];

//truss 3
parameter Real X3=90; //degree between truss
Real k3=A*E/1.25;
Real K3[4,4]; //stiffness matrice
Integer p3a=2;
Integer p3b=4;
Real G3[N,N];

//truss 4
parameter Real X4=90+38.6598; //degree between truss
Real k4=A*E/1.6;
Real K4[4,4]; //stiffness matrice
Integer p4a=1;
Integer p4b=4;
Real G4[N,N];

//truss 5
parameter Real X5=90-38.6598; //degree between truss
Real k5=A*E/1.6;
Real K5[4,4]; //stiffness matrice
Integer p5a=3;
Integer p5b=4;
Real G5[N,N];

/*
for each truss, please ensure pXa is lower then pXb (X represents truss element number)
*/

algorithm

//creating global matrice
K1:=Stiffness_Matrices(X1);
G1:=k1*Local_Global(K1,N,p1a,p1b);

K2:=Stiffness_Matrices(X2);
G2:=k2*Local_Global(K2,N,p2a,p2b);

K3:=Stiffness_Matrices(X3);
G3:=k3*Local_Global(K3,N,p3a,p3b);

K4:=Stiffness_Matrices(X4);
G4:=k4*Local_Global(K4,N,p4a,p4b);

K5:=Stiffness_Matrices(X5);
G5:=k5*Local_Global(K5,N,p5a,p5b);

G:=G1+G2+G3+G4+G5;
Ginitial:=G;

//implementing boundary condition
for i in 1:N loop
 G[2*b1-1,i]:=0;
 G[2*b1,i]:=0;
 G[2*b2-1,i]:=0;
 G[2*b2,i]:=0;
end for;

G[2*b1-1,2*b1-1]:=1;
G[2*b1,2*b1]:=1;
G[2*b2-1,2*b2-1]:=1;
G[2*b2,2*b2]:=1;

//solving displacement
Sol:=Gauss_Jordan(N,G,X);

//solving reaction force
SolMat:=matrix(Sol);
XMat:=matrix(X);
R:=Reaction_Trusses(N,Ginitial,SolMat,XMat);

end Trusses_HW;

Fungsi Panggil

Matrice Transformation

function Stiffness_Matrices
input Real A;
Real Y;
output Real X[4,4];
Real float_error = 10e-10;

final constant Real pi=2*Modelica.Math.asin(1.0);

algorithm

Y:=A/180*pi;
    
X:=[(Modelica.Math.cos(Y))^2,Modelica.Math.cos(Y)*Modelica.Math.sin(Y),-(Modelica.Math.cos(Y))^2,-Modelica.Math.cos(Y)*Modelica.Math.sin(Y);

Modelica.Math.cos(Y)*Modelica.Math.sin(Y),(Modelica.Math.sin(Y))^2,-Modelica.Math.cos(Y)*Modelica.Math.sin(Y),-(Modelica.Math.sin(Y))^2;

-(Modelica.Math.cos(Y))^2,-Modelica.Math.cos(Y)*Modelica.Math.sin(Y),(Modelica.Math.cos(Y))^2,Modelica.Math.cos(Y)*Modelica.Math.sin(Y);

-Modelica.Math.cos(Y)*Modelica.Math.sin(Y),-(Modelica.Math.sin(Y))^2,Modelica.Math.cos(Y)*Modelica.Math.sin(Y),(Modelica.Math.sin(Y))^2];

for i in 1:4 loop
 for j in 1:4 loop
   if abs(X[i,j]) <= float_error then
     X[i,j] := 0;
   end if;
 end for;
end for;

end Stiffness_Matrices;

Global Element Matrice

function Local_Global
input Real Y[4,4];
input Integer B;
input Integer p1;
input Integer p2;
output Real G[B,B];

algorithm

for i in 1:B loop
 for j in 1:B loop
     G[i,j]:=0;
 end for;
end for;

G[2*p1,2*p1]:=Y[2,2];
G[2*p1-1,2*p1-1]:=Y[1,1];
G[2*p1,2*p1-1]:=Y[2,1];
G[2*p1-1,2*p1]:=Y[1,2];

G[2*p2,2*p2]:=Y[4,4];
G[2*p2-1,2*p2-1]:=Y[3,3];
G[2*p2,2*p2-1]:=Y[4,3];
G[2*p2-1,2*p2]:=Y[3,4];

G[2*p2,2*p1]:=Y[4,2];
G[2*p2-1,2*p1-1]:=Y[3,1];
G[2*p2,2*p1-1]:=Y[4,1];
G[2*p2-1,2*p1]:=Y[3,2];

G[2*p1,2*p2]:=Y[2,4];
G[2*p1-1,2*p2-1]:=Y[1,3];
G[2*p1,2*p2-1]:=Y[2,3];
G[2*p1-1,2*p2]:=Y[1,4];

end Local_Global;

Reaction Matrice Equation

function Reaction_Trusses
input Integer N;
input Real A[N,N];
input Real B[N,1];
input Real C[N,1];
Real X[N,1];
output Real Sol[N];
Real float_error = 10e-10;

algorithm
X:=A*B-C;

for i in 1:N loop
 if abs(X[i,1]) <= float_error then
   X[i,1] := 0;
 end if;
end for;

for i in 1:N loop
 Sol[i]:=X[i,1];
end for;

end Reaction_Trusses;

Gauss Jordan

function Gauss_Jordan
input Integer N;
input Real A[N,N];
input Real B[N];
output Real X[N];
Real float_error = 10e-10;

algorithm
X:=Modelica.Math.Matrices.solve(A,B);

for i in 1:N loop
 if abs(X[i]) <= float_error then
   X[i] := 0;
 end if;
end for;

end Gauss_Jordan;

Kuis Membuat Class FLowchart - Metode Numerik/18 November 2020

Kuis Metnum1 Jos.jpegKuis Metnum2 Jos.jpeg

3D Trusses - Metode Numerik/18 November 2020

3D Trusses Flowchart Jos.jpeg
3D Trusses Soal Jos.png
Verifikasi Metnum Jos.jpeg


Master

class Trusses3D_Tes

//define initial variable
parameter Integer Points=4; //Number of Points
parameter Integer Trusses=3; //Number of Trusses
parameter Real Area=0.0015; //Area
parameter Real Elas=70e9; //Elasticity

//define connection
parameter Integer C[Trusses,2]=[1,2;
                                1,3;
                                1,4];
                              
//define coordinates (please put orderly)
parameter Real P[Points,3]=[2,0,0;
                            0,0,1.5;
                            0,0,-1.5;
                            0,1.5,0]; 

//define external force (please put orderly)
parameter Real F[Points*3]={0,-5000,0,
                            0,0,0, 
                            0,0,0, 
                            0,0,0}; 

//define boundary
parameter Integer b[:]={2,3,4};

//solution
Real displacement[N], reaction[N];

protected
parameter Integer N=3*Points;
Integer boundary[3*size(b,1)]=cat(1,(3*b).-2,(3*b).-1,3*b);
Real q1[3], q2[3], g[N,N], G[N,N], G_star[N,N], id[N,N]=identity(N), err=10e-10, cx, cy, cz, L, X[3,3];

algorithm
//Creating Global Matrix
G:=id;
for i in 1:Trusses loop
 for j in 1:3 loop
   q1[j]:=P[C[i,1],j];
   q2[j]:=P[C[i,2],j];
 end for;

       //Solving Matrix
       L:=Modelica.Math.Vectors.length(q2-q1);
       cx:=(q2[1]-q1[1])/L;
       cy:=(q2[2]-q1[2])/L;
       cz:=(q2[3]-q1[3])/L; 
       X:=(Area*Elas/L)*[cx^2,cx*cy,cx*cz;
                         cy*cx,cy^2,cy*cz;
                         cz*cx,cz*cy,cz^2];

       //Transforming to global matrix
       g:=zeros(N,N); 
       for m,n in 1:3 loop
         g[3*(C[i,1]-1)+m,3*(C[i,1]-1)+n]:=X[m,n];
         g[3*(C[i,2]-1)+m,3*(C[i,2]-1)+n]:=X[m,n];
         g[3*(C[i,2]-1)+m,3*(C[i,1]-1)+n]:=-X[m,n];
         g[3*(C[i,1]-1)+m,3*(C[i,2]-1)+n]:=-X[m,n];
       end for;  
 
 G_star:=G+g;
 G:=G_star;
end for;

//Implementing boundary
for i in boundary loop
 for j in 1:N loop
   G[i,j]:=id[i,j];
 end for;
end for;

//Solving displacement
displacement:=Modelica.Math.Matrices.solve(G,F);

//Solving reaction
reaction:=(G_star*displacement)-F;

//Eliminating float error
for i in 1:N loop
 reaction[i]:=if abs(reaction[i])<=err then 0 else reaction[i];
 displacement[i]:=if abs(displacement[i])<=err then 0 else displacement[i];
end for;

end Trusses3D_Tes;

Berikut adalah hasil simulasi yang dilakukan:

DisplacementReaction Forces

Berikut link Youtube yang berisi penjelasan tentang algoritma dan penulisan kode OpenModelica metode diatas:

Optimization - Metode Numerik/16 November 2020

Function to be solved

function Func_Optimization
input Real x;
import Modelica.Math;
output Real y;

algorithm
y:=2*sin(x)-x^2/10;

end Func_Optimization;

Optimization - Golden Section

model Opt_Gold

parameter Real xlo=0;
parameter Real xhi=4; 
parameter Integer N=8; // maximum iteration
parameter Real es=0.0001; // maximum error

Real f1[N], f2[N], x1[N], x2[N], ea[N];
Real xopt,  fx;
protected
Real d, xl, xu, xint, R=(5^(1/2)-1)/2;

algorithm
xl := xlo; 
xu := xhi;
 
for i in 1:N loop
 d:= R*(xu-xl);
 x1[i]:=xl+d;
 x2[i]:=xu-d;
 f1[i]:=Func_Optimization(x1[i]);
 f2[i]:=Func_Optimization(x2[i]);
 xint:=xu-xl;
 
 if f1[i]>f2[i] then
   xl:=x2[i];
   xopt:=x1[i];
   fx:=f1[i];
   else
     xu:=x1[i];
     xopt:=x2[i];
     fx:=f2[i];
 end if;
 
 ea[i]:=(1-R)*abs((xint)/xopt);
 if ea[i]<es then
   break;
 end if;
end for;

end Opt_Gold;

Optimization - Newton

model Opt_Newton

parameter Real g=2.5; //initial guess
parameter Integer N=8; // maximum iteration
parameter Real es=0.0001; // maximum error

Real X[N];
Real xopt, fx, ea[N];

algorithm
X[1]:=g;
ea[1]:=1;

for i in 2:N loop
 X[i]:=X[i-1]-Func_Optimization_Der(X[i-1])/Func_Optimization_Der_Der(X[i-1]);
 xopt:=X[i];
 
 ea[i]:=abs(1-X[i-1]/X[i]);
 if ea[i]<es then
   break;
 end if;
end for;
fx:=Func_Optimization(xopt);

end Opt_Newton;

Optimization - Parabolic

model Opt_Parabolic

parameter Real g1=0; //initial guess
parameter Real g2=1; //initial guess
parameter Real g3=4; //initial guess
parameter Integer N=5; // maximum iteration
parameter Real es=0.0001; // maximum error

Real x1, x2, x3, xopt, xp[N], ea[N];
//Real xl, xm, xu;
Real fx1, fx2, fx3, fx, A[4], A_star[4];
algorithm
x1:=g1;
x2:=g2;
x3:=g3;

for i in 1:N loop
 fx1:=Func_Optimization(x1);
 fx2:=Func_Optimization(x2);
 fx3:=Func_Optimization(x3);
 
 xp[i]:=(fx1*(x2^2-x3^2)+fx2*(x3^2-x1^2)+fx3*(x1^2-x2^2))/(2*fx1*(x2-x3)+2*fx2*(x3-x1)+2*fx3*(x1-x2));
 xopt:=xp[i];
 fx:=Func_Optimization(xp[i]);
 A:={x1,x2,x3,xp[i]};
 A_star:=Modelica.Math.Vectors.sort(A);
  
 if xp[i]>x2 then
   x1:=A_star[2];
   x2:=A_star[3];
   x3:=A_star[4];
   else
     x1:=A_star[1];
     x2:=A_star[2];
     x3:=A_star[3];
 end if;
end for;

ea[1]:=1;
for i in 2:N loop  
 ea[i]:=abs(1-xp[i-1]/xp[i]);
 if ea[i]<es then
   break;
 end if;
end for;

end Opt_Parabolic;


Tugas Besar Metode Numerik (Latihan)- Metode Numerik/23 Desember 2020

Tugas Besar Metode Numerik - Josiah Enrico S

Jawaban UAS 2021 - Josiah Enrico S

Revisi dan Narasi

No 1

Langkah-langkah optimasi

  1. Tentukan tujuan optimasi lalu definisikan semua parameter yang ada.
  2. Menentukan boundary dan constraint pada struktur
  3. Menentukan asumsi dalam permodelan struktur
  4. Melakukan analisa struktur (komputasi)
  5. Setelah mendapat persamaan utama, variasikan nilai objektif yang akan dioptimasi
  6. Gunakan metode optimasi yang sesuai
  7. Setelah mendapatakan nilai optimasi, lakukan pengujian variabel dengan analisa
  8. Jika hasilnya memuaskan, optimasi selesai. Jika hasilnya belum memuaskan, ulangi dari langkah 6

No 2

Tujuan Optimasi:
  • Mengoptimasi kekuatan material struktur dan dinding water tank
Asumsi:
  • Fluida Statis
  • Gravitasi: 9,81 m/s2
  • Fluida: air murni (rho = 1000 kg/m3)
Dalil Fisika:
  • Gaya Hidrostatis F = rho*g*V
  • Stress = E*Strain

No 3

Nomor 3 Ilustrasi struktur truss water tower (3 elemen)
Matriks Persamaan Aljabar
Persamaan aljabar lokal tiap elemen Persamaan aljabar global tiap elemen

Truss A (1 to 2)

            [ 0  0  0  0 ] (1x)
k = (A*E/L) [ 0  1  0  1 ] (1y)
            [ 0  0  0  0 ] (2x)
            [ 0 -1  0 -1 ] (2y)

Truss A (1 to 2)

            [ 0  0  0  0  0  0  0  0 ] (1x)
            [ 0  1  0  1  0  0  0  0 ] (1y)
            [ 0  0  0  0  0  0  0  0 ] (2x)
k = (A*E/L) [ 0 -1  0 -1  0  0  0  0 ] (2y)
            [ 0  0  0  0  0  0  0  0 ] (3x)
            [ 0  0  0  0  0  0  0  0 ] (3y)
            [ 0  0  0  0  0  0  0  0 ] (4x)
            [ 0  0  0  0  0  0  0  0 ] (4y)

Truss B (2 to 3)

            [ 0.5  0.5 -0.5 -0.5 ] (2x)
k = (A*E/L) [ 0.5  0.5 -0.5 -0.5 ] (2y)
            [-0.5 -0.5  0.5  0.5 ] (3x)
            [-0.5 -0.5  0.5  0.5 ] (3y)

Truss B (2 to 3)

            [ 0  0  0    0    0    0    0  0 ] (1x)
            [ 0  0  0    0    0    0    0  0 ] (1y)
            [ 0  0  0.5  0.5 -0.5 -0.5  0  0 ] (2x)
k = (A*E/L) [ 0  0  0.5  0.5 -0.5 -0.5  0  0 ] (2y)
            [ 0  0 -0.5 -0.5  0.5  0.5  0  0 ] (3x)
            [ 0  0 -0.5 -0.5  0.5  0.5  0  0 ] (3y)
            [ 0  0  0    0    0    0    0  0 ] (4x)
            [ 0  0  0    0    0    0    0  0 ] (4y)

Truss C (2 to 4)

            [ 0.5 -0.5 -0.5  0.5 ] (2x)
k = (A*E/L) [-0.5  0.5  0.5 -0.5 ] (2y)
            [-0.5  0.5  0.5 -0.5 ] (4x)
            [ 0.5 -0.5 -0.5  0.5 ] (4y)

Truss C (2 to 4)

            [ 0  0  0    0    0  0  0    0   ] (1x)
            [ 0  0  0    0    0  0  0    0   ] (1y)
            [ 0  0  0.5 -0.5  0  0 -0.5  0.5 ] (2x)
k = (A*E/L) [ 0  0 -0.5  0.5  0  0  0.5 -0.5 ] (2y)
            [ 0  0  0    0    0  0  0    0   ] (3x)
            [ 0  0  0    0    0  0  0    0   ] (3y)
            [ 0  0 -0.5  0.5  0  0  0.5 -0.5 ] (4x)
            [ 0  0  0.5 -0.5  0  0 -0.5  0.5 ] (4y)

No 4

No 5

Fungsi Objektif

displacement= [A*E/L][Uglobal]* deltaL = [F]
Reaction Force= R = [A*E/L][Uglobal]*deltaL - [F]
Stress= Stress = E*(deltaL)/A

Constraint/Boundary

  • Berdasarkan struktur di nomor 3, boundary condition terdapat pada titik 3 dan 4

No 6

Asumsi:

  • Area: 0.16 (Persegi 0.4 x 0.4 meter)
  • Elastisitas: 193 GPa (SS 304)
  • Yield: 215 MPa (SS 304)
  • rho air murni: 1000 kg/m3
  • Gravitasi: 9,81 m/s2
  • Volume air: 200000 gallon / 1136 m3

No 7

Kode Modelica

model UAS_METNUM

//define initial variable
parameter Integer Points=4; //Number of Points
parameter Integer Trusses=3; //Number of Trusses
parameter Real Yield=215e6; //Yield Strength (Pa)
parameter Real Area=0.16;   //Area Block Profile (0.4 * 0.4) (m2)
parameter Real Elas=193e9;     //Elasticity SS 304  (Pa)
parameter Real Volume=1135; //m3 from 200000 gallons
parameter Real rho=1000; //kg/m3

parameter Real W=Volume*rho*9.81;

//define connection
parameter Integer C[Trusses,2]=[1,2;
                                2,3;
                                2,4];
                                                             
//define coordinates (please put orderly)
parameter Real P[Points,3]=[0,37,0;
                            0,8,0;
                           -8,0,0;
                            8,0,0];
                           
//define external force (please put orderly)
parameter Real F[Points*3]={0,-W,0,
                            0,0,0, 
                            0,0,0, 
                            0,0,0}; 

//define boundary
parameter Integer b[:]={3,4};

//solution
Real displacement[N], reaction[N];
Real check[3];

Real stress1[Trusses];
Real safety[Trusses];
Real dis[3];
Real Str[3];

protected
parameter Integer N=3*Points;
Integer boundary[3*size(b,1)]=cat(1,(3*b).-2,(3*b).-1,3*b);
Real q1[3], q2[3], g[N,N], G[N,N], G_star[N,N], id[N,N]=identity(N), cx, cy, cz, L, X[3,3];
Real err=10e-10, ers=10e-4;

algorithm
//Creating Global Matrix
G:=id;
for i in 1:Trusses loop
for j in 1:3 loop
  q1[j]:=P[C[i,1],j];
  q2[j]:=P[C[i,2],j];
end for;
      
   //Solving Matrix
   L:=Modelica.Math.Vectors.length(q2-q1);
   cx:=(q2[1]-q1[1])/L;
   cy:=(q2[2]-q1[2])/L;
   cz:=(q2[3]-q1[3])/L; 
   X:=(Area*Elas/L)*[cx^2,cx*cy,cx*cz;
                     cy*cx,cy^2,cy*cz;
                     cz*cx,cz*cy,cz^2];

   //Transforming to global matrix
   g:=zeros(N,N); 
   for m,n in 1:3 loop
     g[3*(C[i,1]-1)+m,3*(C[i,1]-1)+n]:=X[m,n];
     g[3*(C[i,2]-1)+m,3*(C[i,2]-1)+n]:=X[m,n];
     g[3*(C[i,2]-1)+m,3*(C[i,1]-1)+n]:=-X[m,n];
     g[3*(C[i,1]-1)+m,3*(C[i,2]-1)+n]:=-X[m,n];
   end for;  

G_star:=G+g;
G:=G_star;
end for;

//Implementing boundary
for i in boundary loop
for j in 1:N loop
  G[i,j]:=id[i,j];
end for;
end for;

//Solving displacement
displacement:=Modelica.Math.Matrices.solve(G,F);

//Solving reaction
reaction:=(G_star*displacement)-F;

//Eliminating float error
for i in 1:N loop
reaction[i]:=if abs(reaction[i])<=err then 0 else reaction[i];
displacement[i]:=if abs(displacement[i])<=err then 0 else displacement[i];
end for;

//Checking Force
check[1]:=sum({reaction[i] for i in (1:3:(N-2))})+sum({F[i] for i in (1:3:(N-2))});
check[2]:=sum({reaction[i] for i in (2:3:(N-1))})+sum({F[i] for i in (2:3:(N-1))});
check[3]:=sum({reaction[i] for i in (3:3:N)})+sum({F[i] for i in (3:3:N)});
 
for i in 1:3 loop
 check[i] := if abs(check[i])<=ers then 0 else check[i];
end for;

//Calculating stress in each truss
for i in 1:Trusses loop
for j in 1:3 loop
  q1[j]:=P[C[i,1],j];
  q2[j]:=P[C[i,2],j];
  dis[j]:=abs(displacement[3*(C[i,1]-1)+j]-displacement[3*(C[i,2]-1)+j]);
end for;
      
   //Solving Matrix
   L:=Modelica.Math.Vectors.length(q2-q1);
   cx:=(q2[1]-q1[1])/L;
   cy:=(q2[2]-q1[2])/L;
   cz:=(q2[3]-q1[3])/L; 
   X:=(Elas/L)*[cx^2,cx*cy,cx*cz;
                cy*cx,cy^2,cy*cz;
                cz*cx,cz*cy,cz^2];
   
   Str:=(X*dis);
   stress1[i]:=Modelica.Math.Vectors.length(Str);
end for;

//Safety factor
for i in 1:Trusses loop
 if stress1[i]>0 then
   safety[i]:=Yield/stress1[i];
 else
   safety[i]:=0;
 end if; 
end for;

end UAS_METNUM;