rayleigh ritz method .txt

matlab 求解微分方程

clear
clc
disp('This module Does FEM analysis of a field problem when the Governing differential equation of field is known')
disp('GOVERNING DIFFERENTIAL EQUATION d/dx(a*du/dx)+cu+q=0 ')
disp('NOTE : ENTER VALUES MAKING "a" POSITIVE')
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