Thursday 16 June 2016

Finite Element Method: 1D Elasticity (MATLAB)

Code

%code solves 1D elasticity problems with up to 4 elements (5 nodes)

clc;

n=input('Press one for one, two for two, three for three and four for four elements in the system '); %elements in the system

clc;

if n==1

A1=input('Enter the area of the element in mm^2 '); %area input

E1=input('Enter the Young''s modulus for the element in MPa '); %Young's modulus input

L1=input('Enter the length of the element in mm '); %length input

DT=input('Enter the change in temperature in degree Celsius ');

TA1=input('Enter the coefficient of thermal expansion for the element in per degree Celsius ');

Px1i=input('Enter the pressure acting on the face at node i '); Px1j=input('Enter the pressure acting on the face at node j '); %pressure input

X=input('Enter the body force acting on the part under examination in N/mm^3 '); %body force input

P1i=input('Enter the nodal force on node i '); P1j=input('Enter the nodal force on node j '); %nodal force input

clc;

I=[1 -1;-1 1]; kg=((A1*E1)/L1)*I; %global stiffness matrix calculation

FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; %thermal force vector calculation

FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1]; fp1j=A1*Px1j*FPj.'; %pressure force vector

FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; %body force vector

FN1=[P1i P1j].'; %nodal force vector

fg=ft1+fp1i+fp1j+bf1+FN1; %global force vector calculation

u=input('Press one for u1=0, two for u2=0 '); %boundary conditions

clc;

if u==1

fgu=zeros(1,1); fgu(1,1)=fg(2,1); %global force vector when u1=0

kgu=zeros(1,1); kgu(1,1)=kg(2,2); %global stiffness matrix when u1=0

U=kgu\fgu; %displacement calculation

e1=(U(1,1)-0)/L1; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; %strain due to forces

S1=ef1*E1; %stress calculation

R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction force calculation

uU=['Displacement at node j is ',num2str(U),' mm']; disp(uU); %displacement output

rR=['Reaction force at node i is ',num2str(R1),' N']; disp(rR); %reaction force output

eE=['Strain in the element is ',num2str(ef1)]; disp(eE); %strain output

sS=['Stress in the element is ',num2str(S1),' MPa']; disp(sS); %stress output

end

if u==2

fgu=zeros(1,1); fgu(1,1)=fg(1,1); %global force vector when u1=0

kgu=zeros(1,1); kgu(1,1)=kg(1,1); %global stiffness matrix when u1=0

U=kgu\fgu; %displacement calculation

e1=(0-U(1,1))/L1; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; %strain due to forces

S1=ef1*E1; %stress calculation

R1=(kg(2,1)*U(1,1))-fg(2,1); %reaction force calculation

uU=['Displacement at node i is ',num2str(U),' mm']; disp(uU); %displacement output

rR=['Reaction force at node j is ',num2str(R1),' N']; disp(rR); %reaction force output

eE=['Strain in the element is ',num2str(ef1)]; disp(eE); %strain output

sS=['Stress in the element is ',num2str(S1),' MPa']; disp(sS); %stress output

end

if n==2

A1=input('Enter the area of element number one in mm^2 '); A2=input('Enter the area of element number two in mm^2 '); %area input

E1=input('Enter the Young''s modulus for element number one in MPa '); E2=input('Enter the Young''s modulus for element number two in MPa '); %Young's modulus input

L1=input('Enter the length of element one in mm '); L2=input('Enter the length of element two in mm '); %length input

DT=input('Enter the change in temperature in degree Celsius ');

Px1i=input('Enter the pressure acting on the face at node i for element number one '); Px2j=input('Enter the pressure acting on the face at node j for element number two '); %pressure input

X=input('Enter the body force acting on the part under examination in N/mm^3 '); %body force input

clc;

I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I; ke2=((A2*E2)/L2)*I; %element stiffness matrix calculation

kg=zeros(n+1,n+1); %initialize global stiffness matrix

kg(1,1)=ke1(1,1); kg(1,2)=ke1(1,2); kg(1,3)=0; kg(2,1)=ke1(2,1); kg(2,2)=ke1(1,1)+ke2(2,2); kg(2,3)=ke2(2,1); kg(3,1)=0; kg(3,2)=ke2(2,1); kg(3,3)=ke2(2,2); %global stiffness matrix formulation

FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; ft2=E2*A2*TA2*DT*FT.'; %thermal force vector

FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1]; fp2j=A2*Px2j*FPj.'; %pressure force vector

FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; bf2=(X*A2*L2/2)*FB.'; %body force vector

FN1=[P1i P1j].'; FN2=[P2i P2j].'; %nodal force vector

f1=ft1+fp1i+bf1+FN1; f2=ft2+fp2j+bf2+FN2; %force vector calculation

fg=zeros(3,1); %global force vector initialize

fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1); fg(3,1)=f2(2,1); %global force vector calculation

u=input('Press one for u1=0, two for u3=0, three for u1=u3=0 '); %boundary conditions

clc;

if u==1

fgu=zeros(2,1); fgu(1,1)=fg(2,1); fgu(2,1)=fg(3,1); %global force vector when u1=0

kgu=zeros(2,2); kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); %global stiffness matrix when u1=0

U=kgu\fgu; %displacement calculation

e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2; %stress calculation

R1=(kg(2,1)*U(1,1))-fg(1,1); %reaction force calculation

uU=['Displacement at node two and three is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement output

rR=['Reaction force at element one, node i is ',num2str(R1),' N']; disp(rR); %reaction force output

eE=['Strain in the element one and two is ',num2str(ef1),' and ',num2str(ef2)]; disp(eE); %strain output

sS=['Stress in the element one and two is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output

end

if u==2

fgu=zeros(2,1); fgu(1,1)=fg(1,1); fgu(2,1)=fg(2,1); %global force vector when u3=0

kgu=zeros(2,2); kgu(1,1)=kg(1,1); kgu(1,2)=kg(1,2); kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); %global stiffness matrix when u3=0

U=kgu\fgu; %displacement calculation

R1=(kg(3,2)*U(2,1))-fg(3,1); %reaction force calculation

e1=(U(2,1)-U(1,1))/L1; e2=(0-U(2,1))/L2; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2; %stress calculation

uU=['Displacement at node one and two is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement output

rR=['Reaction force at element two, node j is ',num2str(R1),' N']; disp(rR); %reaction force output

eE=['Strain in the element one and two is ',num2str(ef1),' and ',num2str(ef2)]; disp(eE); %strain output

sS=['Stress in the element one and two is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output

end

if u==3

fgu=zeros(1,1); fgu(1,1)=fg(2,1); %global force vector when u3=0

kgu=zeros(1,1); kgu(1,1)=kg(2,2); %global stiffness matrix when u1=u3=0

U=kgu\fgu; %displacement calculation

R1=(kg(2,1)*U(1,1))-fg(1,1); %reaction force calculation

R3=(kg(3,2)*U(1,1))-fg(3,1); %reaction force calculation

e1=(U(1,1)-0)/L1; e2=(0-U(1,1))/L2; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2;%stress calculation

uU=['Displacement at node two is ',num2str(U(1,1)),' mm']; disp(uU); %displacement output

rR1=['Reaction force at element one, node i is ',num2str(R1),' N']; disp(rR1); %reaction force output

rR3=['Reaction force at element two, node j is ',num2str(R3),' N']; disp(rR3); %reaction force output

eE=['Strain in the element one and two is ',num2str(ef1), ' and ',num2str(ef2)]; disp(eE); %strain output

sS=['Stress in the element one and two is ',num2str(S1),' MPa and ',num2str(S2),' MPa']; disp(sS); %stress output

end

if n==3

A1=input('Enter the area of element number one in mm^2 '); A2=input('Enter the area of element number two in mm^2 '); A3=input('Enter the area of element number three in mm^2 '); %area input

E1=input('Enter the Young''s modulus for element number one in MPa '); E2=input('Enter the Young''s modulus for element number two in MPa '); E3=input('Enter the Young''s modulus for element number three in MPa '); %Young's modulus input

L1=input('Enter the length of element one in mm '); L2=input('Enter the length of element two in mm '); L3=input('Enter the length of element three in mm '); %length input

DT=input('Enter the change in temperature in degree Celsius ');

TA1=input('Enter the coefficient of thermal expansion for element number one in per degree Celsius '); TA2=input('Enter the coefficient of thermal expansion for element number two in per degree Celsius '); TA3=input('Enter the coefficient of thermal expansion for element number three in per degree Celsius '); %thermal expansion coefficient input

Px1i=input('Enter the pressure acting on the face at node i for element number one '); Px3j=input('Enter the pressure acting on the face at node j for element number three '); %pressure input

X=input('Enter the body force acting on the part under examination in N/mm^3 '); %body force input

P1i=input('Enter the nodal force on node i for element number one '); P1j=input('Enter the nodal force on node j for element number one '); P2i=P1j; P2j=input('Enter the nodal force on node j for element number two '); P3i=P2j; P3j=input('Enter the nodal force on node j for element number three '); %nodal force input

clc;

I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I; ke2=((A2*E2)/L2)*I; ke3=((A3*E3)/L3)*I; %element stiffness matrix calculation

kg=zeros(n+1,n+1); %initialize global stiffness matrix

kg(1,1)=ke1(1,1); kg(1,2)=ke1(1,2); kg(2,1)=ke1(2,1); kg(2,2)=ke1(2,2)+ke2(1,1); kg(2,3)=ke2(1,2); kg(3,2)=ke2(2,1); kg(3,3)=ke2(2,2)+ke3(1,1); kg(3,4)=ke3(1,2); kg(4,3)=ke3(2,1); kg(4,4)=ke3(2,2); %global stiffness matrix calculation

FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; ft2=E2*A2*TA2*DT*FT.'; ft3=E3*A3*TA3*DT*FT.'; %thermal force vector

FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1]; fp3j=A3*Px3j*FPj.'; %pressure force vector

FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; bf2=(X*A2*L2/2)*FB.'; bf3=(X*A3*L3/2)*FB.'; %body force vector

FN1=[P1i P1j].'; FN2=[P2i P2j].'; FN3=[P3i P3j].'; %nodal force vector

f1=ft1+fp1i+bf1+FN1; f2=ft2+bf2+FN2; f3=ft3+fp3j+bf3+FN3; %force vector calculation

fg=zeros(4,1); %global force vector initialize

fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1); fg(3,1)=f2(2,1)+f3(1,1); fg(4,1)=f3(2,1); %global force vector formulation

u=input('Press one for u1=0, two for u4=0, three for u1=u4=0 '); %boundary conditions

clc;

if u==1

fgu=zeros(3,1); fgu(1,1)=fg(2,1); fgu(2,1)=fg(3,1); fgu(3,1)=fg(4,1); %global force vector when u1=0

kgu=zeros(3,3); kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4); kgu(3,2)=kg(4,3); kgu(3,3)=kg(4,4); %global stiffness matrix when u1=0

U=kgu\fgu; %displacement calculation

R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction force calculation

e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; %stress calculation

uU=['Displacement at node two, three and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm and ',num2str(U(3,1)),' mm']; disp(uU); %displacement output

rR=['Reaction force at element one, node i is ',num2str(R1),' N']; disp(rR); %reaction force output

eE=['Strain in the element one, two and three is ',num2str(ef1),', ',num2str(ef2),' and ',num2str(ef3)]; disp(eE); %strain output

sS=['Stress in the element one, two and three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa']; disp(sS); %stress output

end

if u==2

fgu=zeros(3,1); fgu(1,1)=fg(1,1); fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); %global force vector when u4=0

kgu=zeros(3,3); kgu(1,1)=kg(1,1); kgu(1,2)=kg(1,2); kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); kgu(2,3)=kg(2,3); kgu(3,2)=kg(3,2); kgu(3,3)=kg(3,3); %global stiffness matrix when u4=0

U=kgu\fgu; %displacement calculation

R1=(kg(4,3)*U(3,1))-fg(4,1); %reaction force calculation

e1=(0-U(3,1))/L1; e2=(U(3,1)-U(2,1))/L2; e3=(U(2,1)-U(1,1))/L3; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; %stress calculation

uU=['Displacement at node one, two and three is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm and ',num2str(U(3,1)),' mm']; disp(uU); %displacement output

rR=['Reaction force at element three, node j is ',num2str(R1),' N']; disp(rR); %reaction force output

eE=['Strain in the element one, two and three is ',num2str(ef1),', ',num2str(ef2),' and ',num2str(ef3)]; disp(eE); %strain output

sS=['Stress in the element one, two and three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa']; disp(sS); %stress output

end

if u==3

fgu=zeros(2,1); fgu(1,1)=fg(2,1); fgu(2,1)=fg(3,1); %global force vector when u1=u4=0

kgu=zeros(2,2); kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); %global stiffness matrix when u4=u1=0

U=kgu\fgu; %displacement calculation

R1=(kg(1,2)*U(1,1))-fg(1,1); R2=(kg(4,3)*U(2,1))-fg(4,1); %reaction force calculation

e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; e3=(0-U(2,1))/L3; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; %stress calculation

uU=['Displacement at node two and three is ',num2str(U(1,1)),' mm and ',num2str(U(2,1)),' mm']; disp(uU); %displacement output

rR=['Reaction force at element one, node i and element three, node j is ',num2str(R1),' N and ',num2str(R2),' N']; disp(rR); %reaction force output

eE=['Strain in the element one, two and three is ',num2str(ef1),' , ',num2str(ef2),' and ',num2str(ef3)]; disp(eE); %strain output

sS=['Stress in the element one, two and three is ',num2str(S1),' MPa, ',num2str(S2),' MPa and ',num2str(S3),' MPa']; disp(sS); %stress output

end

if n==4

A1=input('Enter the area of element number one in mm^2 '); A2=input('Enter the area of element number two in mm^2 '); A3=input('Enter the area of element number three in mm^2 '); A4=input('Enter the area of element number four in mm^2 '); %area input

L1=input('Enter the length of element one in mm '); L2=input('Enter the length of element two in mm '); L3=input('Enter the length of element three in mm '); L4=input('Enter the length of element four in mm '); %length input

DT=input('Enter the change in temperature in degree Celsius ');

Px1i=input('Enter the pressure acting on the face at node i for element number one '); Px4j=input('Enter the pressure acting on the face at node j for element number four '); %pressure input

X=input('Enter the body force acting on the part under examination in N/mm^3 '); %body force input

clc;

I=[1 -1;-1 1]; ke1=((A1*E1)/L1)*I; ke2=((A2*E2)/L2)*I; ke3=((A3*E3)/L3)*I; ke4=((A4*E4)/L4)*I; %element stiffness matrix calculation

kg=zeros(n+1,n+1); %initialize global stiffness matrix

FT=[-1 1]; ft1=E1*A1*TA1*DT*FT.'; ft2=E2*A2*TA2*DT*FT.'; ft3=E3*A3*TA3*DT*FT.'; ft4=E4*A4*TA4*DT*FT.'; %thermal force vector

FPi=[1 0]; fp1i=A1*Px1i*FPi.'; FPj=[0 1]; fp4j=A4*Px4j*FPj.'; %pressure force vector

FB=[1 1]; bf1=(X*A1*L1/2)*FB.'; bf2=(X*A2*L2/2)*FB.'; bf3=(X*A3*L3/2)*FB.'; bf4=(X*A4*L4/2)*FB.'; %body force vector

FN1=[P1i P1j].'; FN2=[P2i P2j].'; FN3=[P3i P3j].'; FN4=[P4i P4j].';%nodal force vector

f1=ft1+fp1i+bf1+FN1; f2=ft2+bf2+FN2; f3=ft3+bf3+FN3; f4=ft4+fp4j+bf4+FN4; %force vector calculation

fg=zeros(5,1); %global force vector initialize

fg(1,1)=f1(1,1); fg(2,1)=f1(2,1)+f2(1,1); fg(3,1)=f2(2,1)+f3(1,1); fg(4,1)=f3(2,1)+f4(1,1); fg(5,1)=f4(2,1); %global force vector formulation

u=input('Press one for u1=0, two for u5=0, three for u1=u5=0 '); %boundary conditions0

clc;

if u==1

fgu=zeros(4,1); fgu(1,1)=fg(2,1); fgu(2,1)=fg(3,1); fgu(3,1)=fg(4,1); fgu(4,1)=fg(5,1); %global force vector when u1=0

kgu=zeros(4,4);

kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4); kgu(3,2)=kg(4,3); kgu(3,3)=kg(4,4); kgu(3,4)=kg(4,5); kgu(4,3)=kg(5,4); kgu(4,4)=kg(5,5);%global stiffness matrix when u1=0

U=kgu\fgu; %displacement calculation

R1=(kg(1,2)*U(1,1))-fg(1,1); %reaction force calculation

e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; e4=(U(4,1)-U(3,1))/L4; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; ef4=e4-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; S4=ef4*E4; %stress calculation

uU=['Displacement at node two, three, four and five is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm, ',num2str(U(3,1)),' mm and ',num2str(U(4,1)),' mm']; disp(uU); %displacement output

rR=['Reaction force at element one, node i is ',num2str(R1),' N']; disp(rR); %reaction force output

eE=['Strain in the element one, two, three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and ',num2str(ef4)]; disp(eE); %strain output

sS=['Stress in the element one, two, three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa and ',num2str(S4),' MPa']; disp(sS); %stress output

end

if u==2

fgu=zeros(4,1); fgu(1,1)=fg(1,1); fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); fgu(4,1)=fg(4,1); %global force vector when u5=0

kgu=zeros(4,4);

kgu(1,1)=kg(1,1); kgu(1,2)=kg(1,2); kgu(2,1)=kg(2,1); kgu(2,2)=kg(2,2); kgu(2,3)=kg(2,3); kgu(3,2)=kg(3,2); kgu(3,3)=kg(3,3); kgu(3,4)=kg(3,4); kgu(4,3)=kg(4,3); kgu(4,4)=kg(4,4); %global stiffness matrix when u5=0

U=kgu\fgu; %displacement calculation

R1=(kg(5,4)*U(4,1))-fg(5,1); %reaction force calculation

e1=(U(2,1)-U(1,1))/L1; e2=(U(3,1)-U(2,1))/L2; e3=(U(4,1)-U(3,1))/L3; e4=(0-U(4,1))/L4; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; ef4=e4-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; S4=ef4*E4; %stress calculation

uU=['Displacement at node one, two, three and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm, ',num2str(U(3,1)),' mm and ',num2str(U(4,1)),' mm']; disp(uU); %displacement output

rR=['Reaction force at element four, node j is ',num2str(R1),' N']; disp(rR); %reaction force output

eE=['Strain in the element one, two, three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and ',num2str(ef4)]; disp(eE); %strain output

sS=['Stress in the element one, two, three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa and ',num2str(S4),' MPa']; disp(sS); %stress output

end

if u==3

fgu=zeros(3,1); fgu(1,1)=fg(1,1); fgu(2,1)=fg(2,1); fgu(3,1)=fg(3,1); %global force vector when u1=u5=0

kgu=zeros(3,3);

kgu(1,1)=kg(2,2); kgu(1,2)=kg(2,3); kgu(2,1)=kg(3,2); kgu(2,2)=kg(3,3); kgu(2,3)=kg(3,4); kgu(3,2)=kg(4,3); kgu(3,3)=kg(4,4); %global stiffness matrix when u5=0

U=kgu\fgu; %displacement calculation

R1=(kg(5,4)*U(3,1))-fg(5,1); R2=(kg(1,2)*U(1,1))-fg(5,1); %reaction force calculation

e1=(U(1,1)-0)/L1; e2=(U(2,1)-U(1,1))/L2; e3=(U(3,1)-U(2,1))/L3; e4=(0-U(3,1))/L4; %strain calculation

e0=TA1*DT; %thermal strain

ef1=e1-e0; ef2=e2-e0; ef3=e3-e0; ef4=e4-e0; %strain due to forces

S1=ef1*E1; S2=ef2*E2; S3=ef3*E3; S4=ef4*E4; %stress calculation

uU=['Displacement at node two, three and four is ',num2str(U(1,1)),' mm, ',num2str(U(2,1)),' mm, ',num2str(U(3,1)),' mm']; disp(uU); %displacement output

rR=['Reaction force at element one, node i and element four, node j is ',num2str(R1),' N and ',num2str(R2),' N']; disp(rR); %reaction force output

eE=['Strain in the element one, two, three and four is ',num2str(ef1),', ',num2str(ef2),', ',num2str(ef3),' and ',num2str(ef4)]; disp(eE); %strain output

sS=['Stress in the element one, two, three and four is ',num2str(S1),' MPa, ',num2str(S2),' MPa ',num2str(S3),' MPa and ',num2str(S4),' MPa']; disp(sS); %stress output

end

Sample calculations for problem with one element

With first boundary condition (u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 1

Enter the area of the element in mm^2 20

Enter the Young's modulus for the element in MPa 200000

Enter the length of the element in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for the element in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i 0

Enter the pressure acting on the face at node j 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i 0

Enter the nodal force on node j 10000

Press one for u1=0, two for u2=0 1

Result

Displacement at node j is 0.26 mm

Reaction force at node i is -10000 N

Strain in the element is 0.0025

Stress in the element is 500 MPa

With second boundary condition (u₂=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 1

Enter the area of the element in mm^2 20

Enter the Young's modulus for the element in MPa 200000

Enter the length of the element in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for the element in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i 0

Enter the pressure acting on the face at node j 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i -10000

Enter the nodal force on node j for element number one 0

Press one for u1=0, two for u2=0 2

Result

Displacement at node i is -0.26 mm

Reaction force at node j is 10000 N

Strain in the element is 0.0025

Stress in the element is 500 MPa

Sample calculations for problem with two elements

Problem 5.1a, Fagan pg. 94

With first boundary condition (u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 2

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 10

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number two 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 10000

Press one for u1=0, two for u3=0, three for u1=u3=0 1

Result

Displacement at node two and three is 0.26 mm and 0.77 mm

Reaction force at element one, node i is -10000 N

Strain in the element one and two is 0.0025 and 0.005

Stress in the element one and two is 500 MPa and 1000 MPa

With second boundary condition (u₃=0)

Problem 5.1a, Fagan pg. 94 mirrored

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 2

Enter the area of element number one in mm^2 10

Enter the area of element number two in mm^2 20

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number two 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one -10000

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 0

Press one for u1=0, two for u3=0, three for u1=u3=0 2

Result

Displacement at node one and two is -0.77 mm and -0.26 mm

Reaction force at element two, node j is 10000 N

Strain in the element one and two is 0.005 and 0.0025

Stress in the element one and two is 1000 MPa and 500 MPa

With third boundary condition (u₃=u₁=0)

Problem 5.1b, Fagan pg. 94

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 2

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 10

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number two 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 0

Press one for u1=0, two for u3=0, three for u1=u3=0 3

Result

Displacement at node two is 0.0033333 mm

Reaction force at element one, node i is 266.6667 N

Reaction force at element two, node j is -266.6667 N

Strain in the element one and two is -6.6667e-05 and -0.00013333

Stress in the element one and two is -13.3333 MPa and -26.6667 MPa

Sample calculations for problem with three elements

Problem 5.1f, Fagan pg. 95

With first boundary condition (u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 3

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 20

Enter the area of element number three in mm^2 10

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the change in temperature in degree Celsius 0

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number three 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 20000

Enter the nodal force on node j for element number two 0

Enter the nodal force on node j for element number three 10000

Press one for u1=0, two for u4=0, three for u1=u4=0 1

Result

Displacement at node two, three and four is 1.25 mm, 1.5 mm and 2 mm

Reaction force at element one, node i is -50000 N

Strain in the element one, two and three is 0.0125, 0.0025 and 0.005

Stress in the element one, two and three is 2500 MPa, 500 MPa and 1000 MPa

With second boundary condition (u₄=0)

Problem 5.1f, Fagan pg. 95 mirrored

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 3

Enter the area of element number one in mm^2 10

Enter the area of element number two in mm^2 20

Enter the area of element number three in mm^2 20

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the change in temperature in degree Celsius 0

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number three 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one -10000

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two -20000

Enter the nodal force on node j for element number three 0

Press one for u1=0, two for u4=0, three for u1=u4=0 2

Result

Displacement at node one, two and three is -2 mm, -1.5 mm and -1.25 mm

Reaction force at element three, node j is 50000 N

Strain in the element one, two and three is 0.0125, 0.0025 and 0.005

Stress in the element one, two and three is 2500 MPa, 500 MPa and 1000 MPa

With third boundary condition (u₄=u₁=0)

Problem 5.1f, Fagan pg. 95 with temperature load only, no nodal forces

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 3

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 20

Enter the area of element number three in mm^2 10

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number three 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 0

Enter the nodal force on node j for element number three 0

Press one for u1=0, two for u4=0, three for u1=u4=0 3

Result

Displacement at node two and three is 0.0025 mm and 0.005 mm

Reaction force at element one, node i and element three, node j is 300 N and -300 N

Strain in the element one, two and three is -7.5e-05 , -7.5e-05 and -0.00015

Stress in the element one, two and three is -15 MPa, -15 MPa and -30 MPa

Sample calculations for problem with four elements

With first boundary condition (u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 4

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 15

Enter the area of element number three in mm^2 10

Enter the area of element number four in mm^2 5

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the Young's modulus for element number four in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the length of element four in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number four in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number four 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 400

Enter the nodal force on node j for element number one 200

Enter the nodal force on node j for element number two 400

Enter the nodal force on node j for element number three 200

Enter the nodal force on node j for element number four 400

Press one for u1=0, two for u5=0, three for u1=u5=0 1

Results

Displacement at node two, three, four and five is 0.06 mm, 0.12333 mm, 0.17333 mm and 0.22333 mm

Reaction force at element one, node i is -2400 N

Strain in the element one, two, three and four is 0.0005, 0.00053333, 0.0004 and 0.0004

Stress in the element one, two, three and four is 100 MPa, 106.6667 MPa 80 MPa and 80 MPa

With second boundary condition (u₅=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 4

Enter the area of element number one in mm^2 5

Enter the area of element number two in mm^2 10

Enter the area of element number three in mm^2 15

Enter the area of element number four in mm^2 20

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the Young's modulus for element number four in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the length of element four in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number four in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number four 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one -400

Enter the nodal force on node j for element number one -200

Enter the nodal force on node j for element number two -400

Enter the nodal force on node j for element number three -200

Enter the nodal force on node j for element number four -400

Press one for u1=0, two for u5=0, three for u1=u5=0 2

Results

Displacement at node one, two, three and four is -0.22333 mm, -0.17333 mm, -0.12333 mm and -0.06 mm

Reaction force at element four, node j is 2400 N

Strain in the element one, two, three and four is 0.0004, 0.0004, 0.00053333 and 0.0005

Stress in the element one, two, three and four is 80 MPa, 80 MPa 106.6667 MPa and 100 MPa

With third boundary condition (u₅=u₁=0)

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 4

Enter the area of element number one in mm^2 20

Enter the area of element number two in mm^2 15

Enter the area of element number three in mm^2 10

Enter the area of element number four in mm^2 5

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the Young's modulus for element number three in MPa 200000

Enter the Young's modulus for element number four in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the length of element three in mm 100

Enter the length of element four in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number three in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number four in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number four 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 0

Enter the nodal force on node j for element number three 0

Enter the nodal force on node j for element number four 0

Press one for u1=0, two for u5=0, three for u1=u5=0 3

Results

Displacement at node two, three and four is -0.0058 mm, -0.0002 mm, 0.0032 mm

Reaction force at element one, node i and element four, node j is -132 N and 132 N

Strain in the element one, two, three and four is -0.000158, -4.4e-05, -6.6e-05 and -0.000132

Stress in the element one, two, three and four is -31.6 MPa, -8.8 MPa -13.2 MPa and -26.4 MPa

Sample Calculations for problems involving variable cross-sectional area

Problem 5.1d, Fagan pg. 95

The command window at various steps

Press one for one, two for two, three for three and four for four elements in the system 2

Enter the area of element number one in mm^2 25

Enter the area of element number two in mm^2 15

Enter the Young's modulus for element number one in MPa 200000

Enter the Young's modulus for element number two in MPa 200000

Enter the length of element one in mm 100

Enter the length of element two in mm 100

Enter the change in temperature in degree Celsius 50

Enter the coefficient of thermal expansion for element number one in per degree Celsius 2/1000000

Enter the coefficient of thermal expansion for element number two in per degree Celsius 2/1000000

Enter the pressure acting on the face at node i for element number one 0

Enter the pressure acting on the face at node j for element number two 0

Enter the body force acting on the part under examination in N/mm^3 0

Enter the nodal force on node i for element number one 0

Enter the nodal force on node j for element number one 0

Enter the nodal force on node j for element number two 10000

Press one for u1=0, two for u3=0, three for u1=u3=0 1

Results

Displacement at node two and three is 0.21 mm and 0.55333 mm

Reaction force at element one, node i is -10000 N

Strain in the element one and two is 0.002 and 0.0033333

Stress in the element one and two is 400 MPa and 666.6667 MPa

Computational Fluid Dynamics: Finite Difference Approximations VS Analytical Differentiation (MATLAB)

Description of the Program

The program takes user input for the function of the variable ‘x’ e.g. x², sin(x), cos(x)/x² etc., the point on which the derivative is evaluated and initial and final grid sizes.

The program then calculates its first derivative analytically.

It then calculates finite difference approximations of the first derivative of order h, with point A unknown and the point A known, using various step sizes, with maximum and minimum step size enter by the user.

It then calculates and plots the error between the exact value and the approximations of the first derivative.

Code along with Guidelines

clc; %clears the command window

%data input and analytical differentiation

syms x

fx=input('Enter the function of ''x'' to be differentiated '); %input a function like sin(x) or x.^4 + x.^3

clc; %clears the command window

fxp=(diff (fx,x)); %calculate the derivative analytically

x=input('Enter the point at which the derivative is to be evaluated '); %input the value at which the derivative to be evaluated

clc; %clears the command window

y=inline(fxp,'x'); %convert a statement in to a function so it can be evaluated at some point

z=y(x); %derivative evaluated at the point

dxmin=input('Enter the smallest value of increment ');

clc;

dxmax=input('Enter the largest value of increment ');

clc;

%finite difference approximation

for dx=dxmin:0.001:dxmax; %loop for different values of grid spacing, dx=h1

a=1.5;

h1=a*dx; %h2 in the problem sheet

h2=(2*a*dx); %h3 in the problem sheet

fc=cos(5+h2)/(5+h2).^2; %value of the function at point c

fb=cos(5+h1)/(5+h1).^2; %of the function at point b

fa=cos(5-dx)/(5-dx).^2; %of the function at point a

fi=cos(x)/(x).^2; %of the function at point A

fxpu=(fc+fb-(2*fa))/((2*dx)+h1+h2); %formulation with A unknown

fxpa=(fc+fb+(2*fa)-(4*fi))/((-2*dx)+h1+h2); %formulation with A known

hold on

plot ((z-fxpa)/z,dx,'.k') %error plot for approximation with A known

plot ((z-fxpu)/z,dx,'.r') %error plot for approximation with A unknown

end

%text editing

title('Plot Between Error VS Grid Spacing') %title of the graph

xlabel('Error from the Exact Value'); %x-axis label

ylabel('Step Size/Grid Spacing') %y-axis label

text(0.6, 0.9,'Red: Error plot for approximation with A unknown') %identification

text(0.6, 0.85,'Black: Error plot for approximation with A known') %identification

Calculation

Sample Problem

Command Window after running the code

At Step One

Enter the function of 'x' to be differentiated cos(x)/x.^2

At Step Two

Enter the point at which the derivative is to be evaluated 5

At Step Three

Enter the smallest value of increment 0.001

At Step Four

Enter the largest value of increment 1

Result

Observations and Recommendations

It was observed that the error in the finite difference approximation when the value of the function at point A is unknown is less as compared to the case where value of the function at point A is known.

In physical/real-world problems, the exact derivatives of the functions are not known, therefore the error values between the exact and the finite difference approximation cannot be calculated. As it is also observed that the error between finite difference approximation and the exact value of the derivative decreases as the step size decreases, it is advisable to use a step size of ~1/1000 to make a compromise between CPU/RAM usage, time taken to solve the problem and error.

Karush–Kuhn–Tucker Conditions (MATLAB)

Code

%KKT with 2 variables, 2 linear constraints with one equality and one

%inequality constraint at once or 1 equality and no inequality constraint at once or no

%equality and one inequality constraint at once

clc;

syms x1 x2 g1 h1 u1 v1 s1

%input of a function like (x1 - 1.5).^2 + (x2 - 1.5).^2 for up to two variables x1, x2

fx=input('Enter the function of ''x'' to be minimized of the form (x - 4).^2 + (y - 6).^2 ');

clc;

C=input('Enter the number of constraints in the problem, e.g. 10 for one equality and no inequality constraint ');

clc;

if C==01

%constraints input

g1=input('Enter inequality constraint in the form x1+x2-10 ');

clc;

%Lagrange function calculation

L=fx+(u1*(g1+s1.^2));

%gradient conditions

dx1=(diff (L,x1)); dx2=(diff (L,x2)); du1=(diff (L,u1)); ds1=(diff (L,s1));

%simultaneous solution of gradient conditions

[sols1,solu1,solx1,solx2]=solve(dx1==0,dx2==0,du1==0,ds1==0);

a=zeros(3,4);

for i=1:3

a(i,1)=solx1(i,1);

a(i,2)=solx2(i,1);

a(i,3)=sols1(i,1);

a(i,4)=solu1(i,1);

end

SET1=a(1,:); SET2=a(2,:); SET3=a(3,:);

%feasibility and non-negativity of Lagrange multiplier check

if SET1(1,3)>0 || SET1(1,4)>0

disp('Set 1 '); disp (' x1 x2 s u'); disp (SET1)

end

if SET2(1,3)>0 || SET2(1,4)>0

disp('Set 2 '); disp (' x1 x2 s u'); disp (SET2)

end

if SET3(1,3)>0 || SET3(1,4)>0

disp('Set 3 '); disp (' x1 x2 s u'); disp (SET3)

end

if C==10

%constraints input

h1=input('Enter equality constraint in the form x1+x2-10 ');

clc;

%Lagrange function calculation

L=fx+(v1*h1);

%gradient conditions

dx1=(diff (L,x1)); dx2=(diff (L,x2)); dv1=(diff (L,v1));

%simultaneous solution of gradient conditions

[solv1,solx1,solx2]=solve(dx1==0,dx2==0,dv1==0);

a=zeros(1,3);

a(1,1)=solx1;

a(1,2)=solx2;

a(1,3)=solv1;

disp('Set 1 '); disp (' x1 x2 v'); disp (a);

end

if C==11

%constraints input

g1=input('Enter inequality constraint in the form x1+x2-10 ');

clc;

h1=input('Enter equality constraint in the form x1+x2-10 ');

clc;

%Lagrange function calculation

L=fx+(v1*h1)+(u1*(g1+s1.^2));

%gradient conditions

dx1=(diff (L,x1)); dx2=(diff (L,x2)); du1=(diff (L,u1)); ds1=(diff (L,s1)); dv1=(diff (L,v1));

%simultaneous solution of gradient conditions

[sols1,solu1,solv1,solx1,solx2]=solve(dx1==0,dx2==0,du1==0,ds1==0,dv1==0);

a=zeros(3,5);

for i=1:3

a(i,1)=solx1(i,1);

a(i,2)=solx2(i,1);

a(i,3)=sols1(i,1);

a(i,4)=solu1(i,1);

a(i,5)=solv1(i,1);

end

SET1=a(1,:); SET2=a(2,:); SET3=a(3,:);

%feasibility and non-negativity of Lagrange multiplier check

if SET1(1,3)>0 || SET1(1,4)>0

disp('Set 1 '); disp (' x1 x2 s u v'); disp (SET1)

end

if SET2(1,3)>0 || SET2(1,4)>0

disp('Set 2 '); disp (' x1 x2 s u v'); disp (SET2)

end

if SET3(1,3)>0 || SET3(1,4)>0

disp('Set 3 '); disp (' x1 x2 s u v'); disp (SET3)

end

Example 4.29 Arora pg. 128

Step One

In this step the user enters the objective function:

Command Window

Enter the function of 'x' to be minimized (x1 - 1.5).^2 + (x2 - 1.5).^2

Step Two

In this step the user specifies the type(s) and quantity of constraints the objective function is subjected to.

01 means no equality and one inequality constraint, 10 means one equality and no inequality constraint, 11 means one equality and one inequality constraint.

Command Window

Enter the number of constraints in the problem, e.g. 10 for one equality and no inequality constraint 01

Step Three

In this step the user enter the constraint(s).

Command Window

Enter inequality constraint in the form x1+x2-10 x1 + x2 – 2

Step Four

Output of the program; feasible set with candidate minimum point(s).

Command Window

Set 1

x1 x2 s u

1 1 0 1

Example 4.27 Arora pg. 122

Step One

Command Window

Enter the function of 'x' to be minimized (x1 - 1.5).^2 + (x2 - 1.5).^2

Step Two

Command Window

Enter the number of constraints in the problem, e.g. 10 for one equality and no inequality constraint 10

Step Three

Command Window

Enter equality constraint in the form x1+x2-10 x1 + x2 – 2

Step Four

Command Window

Set 1

x1 x2 v

1 1 1

Modifications

Extension to five or even more variables and more constraints can be made via simple modifications to the code. The structure of the new code will remain similar to the code mention above, though it will have many more lines of instructions.

My Blogger Profile

Thursday 16 June 2016

With first boundary condition (u1=0)

With second boundary condition (u2=0)

With first boundary condition (u1=0)

With second boundary condition (u3=0)

With third boundary condition (u3=u1=0)

With first boundary condition (u1=0)

With second boundary condition (u4=0)

With third boundary condition (u4=u1=0)

With first boundary condition (u1=0)

With second boundary condition (u5=0)

With third boundary condition (u5=u1=0)

Example 4.29 Arora pg. 128

Step One

Step Two

Step Three

Step Four

Example 4.27 Arora pg. 122

Step One

Step Two

Step Three

Step Four

Modifications

With first boundary condition (u₁=0)

With second boundary condition (u₂=0)

With first boundary condition (u₁=0)

With second boundary condition (u₃=0)

With third boundary condition (u₃=u₁=0)

With first boundary condition (u₁=0)

With second boundary condition (u₄=0)

With third boundary condition (u₄=u₁=0)

With first boundary condition (u₁=0)

With second boundary condition (u₅=0)

With third boundary condition (u₅=u₁=0)