瑞利里兹法

clear；

clc；

syms x

a=input('Enter value of "a" : ');

c=input('Enter value of "c" : ');

q=input('Enter value of "q" : ');

syms le xa xb Q1e Q2e

w1 = 1-(x-xa)/le;

w2 = (x-xa)/le;

syms u1 u2

u = (1-(x-xa)/le)*u1+(x-xa)/le*u2;

f = diff(w1,x)*a*diff(u,x);

tdwdu1 = int(f,x,xa,xa+le);

f = w1*q;

twq1 = int(f,x,xa,xa+le);

f = c*u*w1;

tcu1 = int(f,x,xa,xa+le);

f = diff(w2,x)*a*diff(u,x);

tdwdu2 = int(f,x,xa,xa+le);

f = w2*q;

twq2 = int(f,x,xa,xa+le);

f = c*u*w2;

tcu2 = int(f,x,xa,xa+le);

a1=tdwdu1/(u1-u2);

% a1 = a11 of EME and - a1 = a12 of EME

a2=tdwdu2/(u1-u2);

% a2 = a21 of EME and - a2 = a22 of EME

m1=[simplify(a1) simplify(-a1);simplify(a2) simplify(-a2)];

%m1 is matrix representing primary variable.

m3=[simplify(twq1);simplify(twq2)];

% m3 is matrix which is on rhs side along with Q; internal source

r11=simplify(tcu1);

r21=collect(r11,u1);

r11=collect(r21,u2);

r12=simplify(tcu2);

r22=collect(r12,u1);

r12=collect(r22,u2);

m4=zeros(2,2);

m5=[Q1e;Q2e];

p=subs(r11,u2,0);

z11=subs(p,u1,1);

p=subs(r11,u1,0);

z12=subs(p,u2,1);

p=subs(r12,u1,0);

z22=subs(p,u2,1);

p=subs(r12,u2,0);

z21=subs(p,u1,1);

m4=[z11 z12;z21 z22];

% m4 is matrix which also contains primary variables and comes along side m1

clc

disp('INFORMATION ABOUT DISCRETIZATION')

N=input('Enter number of elements: ');

n=input('Enter number of nodes: ');

k=zeros(n);

q=zeros(n,1);

disp('Give information about the elements ')

for(i=1:N)

l=input(['\nleft node for element ',num2str(i),': ']);

r=input(['\nright node for element ',num2str(i),': ']);

Xa=input(['\nEnter x coordinate of left node for element ',num2str(i),': ']);

Le=input(['\nEnter length of element ',num2str(i),': ']);

var1=m1(1,1);

var2=subs(var1,xa,Xa);

var1=subs(var2,le,Le);

k(l,l)=k(l,l)+var1;

var1=m1(1,2);

var2=subs(var1,xa,Xa);

var1=subs(var2,le,Le);

k(l,r)=k(l,r)+var1;

var1=m1(2,1);

var2=subs(var1,xa,Xa);

var1=subs(var2,le,Le);

k(r,l)=k(r,l)+var1;

var1=m1(2,2);

var2=subs(var1,xa,Xa);

var1=subs(var2,le,Le);

k(r,r)=k(r,r)+var1;

var1=subs(z11,xa,Xa);

var2=subs(var1,le,Le);

k(l,l)=k(l,l)-var2;

var1=subs(z12,xa,Xa);

var2=subs(var1,le,Le);

k(l,r)=k(l,r)-var2;

var1=subs(z21,xa,Xa);

var2=subs(var1,le,Le);

k(r,l)=k(r,l)-var2;

var1=subs(z22,xa,Xa);

var2=subs(var1,le,Le);

k(r,r)=k(r,r)-var2;

var1=m3(1,1);

var2=subs(var1,xa,Xa);

var1=subs(var2,le,Le);

q(l,1)=q(l,1)+var1;

var1=m3(2,1);

var2=subs(var1,xa,Xa);

var1=subs(var2,le,Le);

q(r,1)=q(r,1)+var1;

end

k

k1=k;

clc

disp('INFORMATION ABOUT SOURCE TERMS\nPress 0 for unknown terms')

disp('Boundary conditions like Du(0),Du(4),etc to be entered here')

for(i=1:n)

f(i,1)=input(['enter value of external source at node ',num2str(i),': ']);

end

disp('global source term matrix ');

f2=f+q

f1=f2;

clc

disp('INFORMATION ABOUT BOUNDARY CONDITIONS')

disp('Boundary conditions like u(0),u(4),etc to be entered here')

u1=zeros(n,1);

p=n;

r=1;

i=n;

c=zeros(1,n);

while(i>=1)

d=input(['Is value of Primary Variable at node ',num2str(i),' known? press 1 for yes 0 for no ']);

if(d==1)

c(1,r)=i;

u1(i,1)=input('enter its value: ');

e(1,r)=u1(i,1);

j=p;

while(j>=1)

if(j~=i)

f1(j,1)=f1(j,1)-k1(j,i)*u1(i,1);

end

j=j-1;

end

k1(i,:)=[];

k1(:,i)=[];

f1(i,:)=[];

r=r+1;

p=p-1;

end

i=i-1;

end

u1=inv(k1)*f1;

clear i;

u=i*ones(n,1);

for(j=1:r-1)

pos=c(1,j);

u(pos,1)=e(1,j);

end

j=1;

for(s=1:n)

if(u(s,1)==i)

u(s,1)=str2num(char(u1(j,1)));

j=j+1;

end

end

clc

disp('THE SOLUTION IS DONE !!! ')

disp('Global Primary Variable matrix')

u

figure (1)

plot(u)

axis tight

xlabel('Nodes')

ylabel('Primary Variable')

title('Variation of primary variables with nodes')

grid

disp('Value of Unknown Secondary Variables')

for(l=1:r-1)

sum=0;

p=c(1,l);

for(j=1:n)

sum=sum+k(p,j)*u(j,1);

end

disp(['Q ',num2str(p),' = ',num2str(sum)])

end