Hola, es verdad tienes razón, el código es el siguiente,
function FEM_code__Coupled()
clear all;close all; clc
n = 5;
nn = n+1;
lgth = 1;
he = lgth/n;
x = (0:he:lgth);
AC = 0.00005;
U= zeros(nn,1);
V=zeros(nn,1);
U(1)=0 ; V(1)=1;
U(nn)=1; V(nn)=exp(1);
c = 1.0;
count = 0;
tic
while(c>0)
[U1,V1] = assembly(U,V,n,he);
c =0.0;
for i=1:nn
if (abs(U(i)-U1(i))>AC)
c= c+1;
break;
end
if (abs(V(i)-V1(i))>AC)
c= c+1;
break ;
end
end
U =U1;
V = V1;
count = count +1;
end
disp('hence solution:');
% saida de variaveis primaria e sequandarias.
diff = abs(U-(x.^2)');
diff1 = abs(V-exp(x)');
fprintf(' No de elementos = %d\n',n);
disp(' x u(FEM) u(Exact) error v(FEM) v(Exact) Error')
disp([x',U,(x.^2)',diff,V,exp(x)',diff1])
fprintf('No de iteracoes=%d\n ',count);
%grafica variaveis primarias
figure(1)
plot(x,U,'--rs','Linewidth',2)
xlabel('x');
ylabel('u(x)');
title('Solucao grafica para o pvc')
figure(2)
plot(x,V,'--rs','Linewidth',2)
xlabel('x')
ylabel('v(x)')
title('Solucao grafica para o PVC')
toc
end
%derivacao de elemntos da matriz e emsamble
function [U1,V1] = assembly(U,V,n,he)
nn=n+1;
L11=zeros(nn,nn); L12=zeros(nn,nn); R1=zeros(nn,1);
L21=zeros(nn,nn); L22=zeros(nn,nn); R2=zeros(nn,1);
lmm=[];
for i=1:n
lmm=[lmm;[i,i+1]];
end
for i=1:n
lm=lmm(i,:);
%geracao de elemnto da matriz k11 e RMS matriz (f1)
% Integracao usando quadratura gaussiana
[k111,k112,k113,k311,k312,k313,k314,f11,f12] = linearelement(he,i);
% Para primeira equacao
k11 = -k111;
k12 = k313*V(lm(1))+k314*V(lm(2))+(k311*U(lm(1))+k312*U(lm(2)));
% para sequnda equcao
k21 = k112;
k22 = -k111-k113;
f1 = f11;
f2 = f12;
% matriz de conetividade assembly
L11(lm,lm) = L11(lm,lm)+k11;
L12(lm,lm) = L12(lm,lm)+k12;
L21(lm,lm) = L21(lm,lm)+k21;
L22(lm,lm) = L22(lm,lm)+k22;
R1(lm) = R1(lm)+f1;
R2(lm) = R2(lm) +f2;
end
% imposicao das condicoes de fronteira
L11(1,:)=0.0; L12(1,:)=0.0;
L11(nn,:)=0.0; L12(nn,:)=0.0;
L21(1,:)=0.0 ; L22(1,:)=0.0;
L21(n,:)=0.0 ; L22(n,:)=0.0;
L11(1,1)=1; L11(nn,nn)=1;
L22(1,1)=1; L22(nn,nn)=1;
R1(1,1)=U(1); R1(nn,1)=U(nn);
R2(1,1)=V(1); R2(nn,1)=V(nn);
% solucao do sistema de equacoes F1
K= [L11,L12;L21,L22];
R=[R1;R2];
d=K\R;
U1= d(1:nn);
V1= d(nn+1:2*nn);
end
function [k111,k112,k113,k311,k312,k313,k314,f11,f12] = linearelement(he,i);
gL = [-sqrt(3/5);0;sqrt(3/5)];
gW = [5/9,8/9,5/9];
k111=zeros(2); k112=zeros(2); k113=zeros(2); k221=zeros(2); k222=zeros(2);k311=zeros(2); k312=zeros(2);
k313=zeros(2); k314=zeros(2); f11=zeros(2,1); f12=zeros(2,1);
for k=1:length(gW)
s = gL(k); w = gW(k);
n = [(1-s)/2, (1+s)/2];
dns = [-1/2, 1/2];
coord = [(i-1)*he;i*he];
x = n*coord(:,1);
J=he/2;
dx=(1/J)*dns;
k111=k111+J*w*dx'*dx;
k112=k112+J*w*n'*dx;
k113=k113+J*w*n'*n;
k221=k221+J*w*n'*dx*n(1);
k222=k222+J*w*n'*dx*n(2);
% Cuando digito los siguientes parametros me sale matriz no singular
% % falta definir
% k311=k311+J*w*n'*dx*n(1);
% k312=k312+J*w*n'*dx*n(2);
% k313=k313+J*w*n'*dx*n(1);
% k314=k314+J*w*n'*dx*n(2);
%
% f11=f11+J*w*n'*(x.^2)';
% f12=f12+J*w*n'*exp(x)';
%
end
end
