problem in finding inverse
3 views (last 30 days)
Show older comments
%function[el5,f0]=data_cholesky(sigma,nel,l)
clc
clear all
sigma=0.15;
l=300;
nel=200; % number of elements
nnel=2; % number of nodes per element
ndof=1; % number of dofs per node
nnode=(nnel-1)*nel+1; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
index=zeros(nnel*ndof,1);
mm=zeros(sdof, sdof);
bcdof(1)=1; % first dof (deflection at left end) is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=2; % 12th dof (slope at the symmetric end) is constrained
bcval(2)=0; % whose described value is 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = zeros(nel+1, 1); % Replace with your initial guess for x
i=1; %boundary condition
T(i) = 300;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
keq=zeros(nel+1,nel+1,l);
keq1=zeros(nel+1,1,l);
keq2=zeros(nel+1,nel,l);
keq3=zeros(nel,nel,l);
A1 =zeros(nel+1,1,l);
A2 =zeros(nel,1,l);
el40 = zeros(nel, l); % Initialize el4 as a 2D matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=nel;
i=1:a;
b0=11;
mean0=zeros(3,b0);
vari0=zeros(3,b0);
kk(:,:,:)=zeros(sdof, sdof,l);
ss(:,:)=zeros(sdof, sdof);
m=zeros(2,2,l);
s=zeros(2,1,l);
k0 = zeros(2,1,l);
el0=zeros(nel,b0,l);
%for b0=1:b0
b=b0-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ka = 200 ; % conductivity
tleng=50*10^-3; % total length
d =1*10^-3;
le=tleng/(nel);
k0 = (ka/le)*[1 -1;-1 1];%
P =3.14*d ; % perimeter;
h = 20; % convective heat transfer
Ac = 3.14*d^2/4; % cross area crosssection
q0 = 0;
m = (P * h * le) / (6 * Ac) * [2 1; 1 2];%
Tinf=30;
s = (q0 + (P*h/Ac)*Tinf)*[le*0.5;le*0.5]; %
s=diag(s);
s = repmat(s, [1, 1, l]);
m = repmat(m, [1, 1, l]);
%k0=repmat(k0, [1, 1, l]);
Q=zeros(nel+1,1);
Q(1)=1;
Q(nel+1)=0;
for iel=1:nel % loop for the total number of elements
edof = nnel*ndof;
start = (iel-1)*(nnel-1)*ndof;
%
for i=1:edof
index(i)=start+i;
end
edof = length (index);
for i=1:edof
ii=index(i) ;
for j=1:edof
jj=index(j) ;
mm(ii,jj)=mm(ii,jj)+m(i,j,l);
ss(ii,jj)=ss(ii,jj)+s(i,j,l);
end
end
end
for l=1:l
d1=tleng/nel;
d=1;
i=1:nel;
j=1:nel;
ee=zeros(nel,nel);
R=zeros(nel,nel);
e(i)=tleng/2*nel+(tleng/nel)*(i-1);
for i=1:nel
for j=1:nel
ee(i,j)=e(i)-e(j);
R(i,j)=sigma^2*exp(-abs(ee(i,j)/d));
end
end
C1=chol(R,'lower');
Z1=randn(nel,1);
alpha=C1*(Z1);
E0=1;
el4= ((1+alpha));
el40(:, l)= el4;
el=el4(1:nel);
el2= 1 + zeros( 1, 100 );
el3=[0.90991 0.98608 0.90314 1.1458 1.0554 0.91364 1.3918 0.61164 0.74829 1.369];
el5=E0*ones(1,20);
for iel=1:nel % loop for the total number of elements
edof = nnel*ndof;
start = (iel-1)*(nnel-1)*ndof;
%
for i=1:edof
index(i)=start+i;
end
c0(iel)=el4(iel);
k(:,:,iel)=c0(iel)*k0;
edof = length (index);
for i=1:edof
ii=index(i) ;
for j=1:edof
jj=index(j) ;
kk(ii,jj,l)=kk(ii,jj,l)+k(i,j,iel);
end
end
end
sss=diag(ss);
keq=kk+mm;
keq1=T(1)*keq(:,1);
keq2 =keq(:,2:end,:);
keq3=keq2(2:end,:,:);
A1(:,:,l)=sss+Q-keq1;
A2(:,:,l)=A1(2:end,:,l);
%A3(:,:,l)=inv(keq3(:,:,l));
%temp(:,:,l)=A3(:,:,l)*A2(:,:,l);
temp(:,:,l)=keq3(:,:,l)\A2(:,:,l);
end
% Define the length and number of elements
temp1=squeeze(temp);
T0 =300 * ones(1,l);
temp1=[T0;temp1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%temp1([2:4],:)=[]
x= linspace(0,tleng,nel+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
x=x';
for i=1:l
plot(x,temp1(:,i))
hold on
end
mean_temp =mean(temp1,2) ;
mean_temp(1:3)
the variation of temp vs x should be continously dectreasing. but when i run it several times .sometimes temp increses with x please help
1 Comment
Image Analyst
on 20 Dec 2023
I have no idea what this weakly commented code does. Try adding more comments to explain the various things. Your subject says "problem in finding inverse" so I guess you want to find the inverse of a matrix. About all I can say is that if you want the inverse of a matrix, try inv
help inv
Answers (0)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!