problem in finding inverse
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
%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 Kommentar
Image Analyst
am 20 Dez. 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
Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!