problem in finding inverse
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
%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 commentaire
Image Analyst
le 20 Déc 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
Réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!