what is finite elemnt method code for Sytem of ODE

2 visualizaciones (últimos 30 días)
Tadese
Tadese el 22 de Oct. de 2023
Movida: John D'Errico el 5 de Nov. de 2023
what is the math lab code of the following by finite element method please please
f"'+(f+g)f''-Kf^5-f'^2-Mf'+lambda*(theta+N*phi)=0------ (1)
g'''+(f+g)g''-kg^5+M*g'=0 ----- (2)
theta''+Pr*(f+g)theta'+Pr*Nb*theta'phi'+Pr*Nt*theta'^2-deltatPr(f+g)(f'+g')theta')+(theta'')(f+g)^2=0------------------------(2)
phi"+(Nt/Nb)*theta"+Sc(f+g)-deltac*Sc[(f+g)^2phi"+(f+g)(f'+g')ph']= 0------(3)
'f','g','theta','phi' are the dependent variables and η is the independent variable
And the boundary conditions are as follows:
f′(η) = 1 + αf″(η)) + βf‴(η) + cf⁗(η) + δf(v)(η),
g(η) =c + φg″(η) + χg‴(η) + ψg⁗(η) + ωg(v)(η),
θ′(η) = - βi1(1 - θ(η)),
ϕ′(η) =- βi2(1 - ϕ(η)), at  η = 0, f′(∞) ⟶ 0, f″(∞) ⟶ 0, f‴(∞) ⟶ 0, g′(∞) ⟶ 0, g″(∞) ⟶ 0, g‴(∞) ⟶ 0, θ(∞) ⟶ 0, ϕ(∞) ⟶ 0
So solving these equations numerically I have used finite element method.

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

Torsten
Torsten el 22 de Oct. de 2023
Movida: Torsten el 22 de Oct. de 2023
First correct your boundary conditions.
You need 3 + 3 + 2 + 2 = 10, and they can only have derivatives up to (order of equations - 1).
Thus f''', f'''', g''', g'''' are not allowed.
After having done this, I suggest you start setting up your problem in bvp4c first to have a reference solution.
  4 comentarios
Mostrar 2 comentarios más antiguos
Tadese
Tadese el 5 de Nov. de 2023
Movida: John D'Errico el 5 de Nov. de 2023
what is the math lab code of the following by finite element method please please
f"'+(f+g)f''-Kf'''''-f'^2-Mf'+lambda*(theta+N*phi)=0------ (1)
g'''+(f+g)g''-kg'''''+M*g'=0 ----- (2)
theta''+Pr*(f+g)theta'+Pr*Nb*theta'phi'+Pr*Nt*theta'^2-deltatPr(f+g)(f'+g')theta')+(theta'')(f+g)^2=0------------------------(2)
phi"+(Nt/Nb)*theta"+Sc(f+g)-deltac*Sc[(f+g)^2phi"+(f+g)(f'+g')ph']= 0------(3)
'f','g','theta','phi' are the dependent variables and η is the independent variable
And the boundary conditions are as follows:
f(0)=g(0)=0,
f′(η) = 1 + αf″(η)) + βf‴(η) + cf⁗(η) + δf'''''(η),
g(η) =c + φg″(η) + χg‴(η) + ψg⁗(η) + ωg'''''(η),
Nb*phi(0)+Nt*theta'(0)=0,
θ′(η) = - βi1(1 - θ(η)), at  η = 0, f′(∞) ⟶ 0, f″(∞) ⟶ 0, f‴(∞) ⟶ 0, g′(∞) ⟶ 0, g″(∞) ⟶ 0, g‴(∞) ⟶ 0, θ(∞) ⟶ 0, ϕ(∞) ⟶ 0
So solving these equations numerically I have used finite element method.
By using bvp4c I can but my dissertation by using Finite element method
so, please please any one how know this math lab code help me. Thank you
This is my math lab code can you help me please
function couple130m
format long
clear all
clc
n=200;
nn=n*2+1;
lgth=6 ; %length of the domain
he=lgth/n; % increment
x = [0.0:he/2:lgth]; % nodal points in the domian
AC = 0.000001;% FOR ACCURACY(tolerence)
K=0.9;
M=1.2;
delta1=0.2;
delta2=0.2;
lamda=0.5;
N=0.5; Nt=0.1;
Pr=0.72;
Sc=0.8;
gamma=0.1;
alpha=-1;
delta=1;
beta1=0.1;
beta2=0.3;
Bi=0.5;
Nb=0.1;
for k=[0.2 0.4 0.6 0.8]
F = sparse(nn+1,1);
G = sparse(nn+2,1);
E = sparse(nn+1,1);
Q = sparse(nn+2,1);
T = sparse(nn+2,1);
H = sparse(nn+2,1);
c = 1.0;
count=0;
tic
while(c>0)
[F1,G1,E1,Q1,T1,H1] = assembly(F,G,E,Q,T,H,n,nn,he,Pr,M,Sc,gamma,...
lamda,delta,Bi,alpha,beta1,beta2,delta1,delta2,Nt,Nb,N,K);
c = 0.0;
for i=1:nn
if(abs(F(i)-F1(i))>AC)
c = c+1;
break;
end
if(abs(G(i)-G1(i))>AC)
c = c+1;
break;
end
if(abs(E(i)-E1(i))>AC)
c = c+1;
break;
end
if(abs(Q(i)-Q1(i))>AC)
c = c+1;
break;
end
if(abs(T(i)-T1(i))>AC)
c = c+1;
break;
end
if(abs(H(i)-H1(i))>AC)
c = c+1;
break;
end
end
F = F1;
G = G1;
E = E1;
Q = Q1;
T = T1;
H = H1;
count=count+1;
end
disp('Hence Solution=:');
x = [0.0:he/2:lgth];
figure (1)
% plot(x,G(1:nn),'--b')
% xlabel('\eta')
% ylabel('g(\eta)')
% hold on
if Nb==0.02
plot(x,G(1:nn),'b','linewidth',0.8)
elseif Nb ==0.04
plot(x,G(1:nn),'--r','linewidth',0.8)
elseif Nb==0.06
plot(x,G(1:nn),'g','linewidth',0.8)
elseif Nb==0.08
plot(x,G(1:nn),'k','linewidth',0.8)
else
plot(x,G(1:nn),'k','linewidth',0.8)
end
xlabel('\eta')
ylabel('f''(\eta)')
hold on;
legend('Nb= 0.02','Nb= 0.04','Nb= 0.06','Nb=0.08')
% figure (2)
% if Nb==0.02
% plot(x,Q(1:nn),'b','linewidth',0.8)
% elseif Nb==0.04
% plot(x,Q(1:nn),'--r','linewidth',0.8)
% elseif Nb==0.06
% plot(x,Q(1:nn),'g','linewidth',0.8)
% elseif Nb==0.08
% plot(x,G(1:nn),'k','linewidth',0.8)
% else
% plot(x,Q(1:nn),'k','linewidth',0.8)
%
% end
% xlabel('\eta')
% ylabel('g''(\eta)')
% hold on;
% legend('Nb= 0.02','Nb= 0.04','Nb= 0.06','Nb=0.08')
% figure (3)
% if Nb==0.02
% plot(x,T(1:nn),'b','linewidth',0.8)
% elseif Nb==0.04
% plot(x,T(1:nn),'r','linewidth',0.8)
% elseif Nb==0.06
% plot(x,T(1:nn),'g','linewidth',0.8)
% elseif Nb==0.08
% plot(x,T(1:nn),'k','linewidth',0.8)
% else
% plot(x,T(1:nn),'k','linewidth',0.8)
% end
% xlabel('\eta')
% ylabel('\theta(\eta)')
% hold on;
% legend('Nb= 0.02','Nb= 0.04','Nb= 0.06','Nb=0.08')
% figure (4)
% if Nb==0.02
% plot(x,H(1:nn),'b','linewidth',0.8)
% elseif Nb==0.04
% plot(x,H(1:nn),'r','linewidth',0.8)
% elseif Nb==0.06
% plot(x,H(1:nn),'g','linewidth',0.8)
% elseif Nb==0.08
% plot(x,H(1:nn),'k','linewidth',0.8)
% else
% plot(x,H(1:nn),'k','linewidth',0.8)
% end
% xlabel('\eta')
% ylabel('\phi(\eta)')
% hold on;
% legend('Nb= 0.02','Nb= 0.04','Nb= 0.06','Nb=0.08')
end
Skin_friction =-E(nn+1)
%Skin_friction =-E(nn+1)
Heat_Transfer =-T(nn+1)
sher_wood =-Q(nn+1)
fprintf('%10.6f\t %10.6f\t \n',full(-E(nn+1)),full(-T(nn+1)),full(-H(nn+1)));
end
function [F1,G1,E1,Q1,T1,H1] = assembly(F,G,E,Q,T,H,n,nn,he,Pr,M,Sc,gamma,...
lamda,delta,Bi,alpha,beta1,beta2,delta1,delta2,Nt,Nb,N,K)
% INITIALIZE MATRICES
l1 = sparse(nn,nn);l2 = sparse(nn,nn);l3 = sparse(nn,nn);l4 = sparse(nn,nn);
l5 = sparse(nn,nn);l6 = sparse(nn,nn);l7 = sparse(nn,nn);l8 = sparse(nn,nn);
l9 = sparse(nn,nn);l10 = sparse(nn,nn);l11 = sparse(nn,nn);l12 = sparse(nn,nn);
l13 = sparse(nn,nn);l14 = sparse(nn,nn);l15 = sparse(nn,nn);l16 = sparse(nn,nn);
l17 = sparse(nn,nn);l18 = sparse(nn,nn);l19 = sparse(nn,nn);l20 = sparse(nn,nn);
l21 = sparse(nn,nn);l22 = sparse(nn,nn);l23 = sparse(nn,nn);l24 = sparse(nn,nn);l25 = sparse(nn,nn);
l26 = sparse(nn,nn);l27 = sparse(nn,nn);l28 = sparse(nn,nn);l29 = sparse(nn,nn);l30 = sparse(nn,nn);
l31 = sparse(nn,nn);l32 = sparse(nn,nn);l33 = sparse(nn,nn);l34 = sparse(nn,nn);l35 = sparse(nn,nn);l36 = sparse(nn,nn);
R1 = sparse(nn,1);R2 = sparse(nn,1);R3 = sparse(nn,1);R4 = sparse(nn,1);R5 = sparse(nn,1);R6 = sparse(nn,1);
%step3: Assembly of element equation i.e.[k]{y}={q}
q11=0.0;
q12=0.0;
q13=0.0;
q21=0.0;
q22=0.0;
q23=0.0;
q31=0.0;
q32=0.0;
q33=0.0;
q41=0.0;
q42=0.0;
q43=0.0;
q51=0.0;
q52=0.0;
q53=0.0;
q61=0.0;
q62=0.0;
q63=0.0;
f1=[q11;q12;q13];
f2=[q21;q22;q23];
f3=[q31;q32;q33];
f4=[q41;q42;q43];
f5=[q51;q52;q53];
f6=[q61;q62;q63];
p=[];
for i=1:n
j=(i-1)*2+1;
p=[p;[j,j+1,j+2]];%Connectivity matrix
end
lmm=[];
for i=1:n
lmm =[lmm;[p(i,1),p(i,2),p(i,3)]];
end
for i=1:n
lm=lmm(i,:);
[kv11,kv1,kv,kv1pa,kv1pb,kv1pc,kvpa,kvpb,kvpc,...
kv1pp11,kv1pp12,kv1pp13,kv1pp21,kv1pp22,kv1pp23,kv1pp31,...
kv1pp32,kv1pp33,kv11pp11,kv11pp12,kv11pp13,kv11pp21,kv11pp22,...
kv11pp23,kv11pp31,kv11pp32,kv11pp33,...
kv1p1a,kv1p1b,kv1p1c,kt,knn] = quadraticelement(he,i);
%F
k11=kv1;
k12=sparse(3,3);
k13=-kv;
k14=sparse(3,3);
k15=sparse(3,3);
k16=sparse(3,3);
%S
k21=sparse(3,3);
k22=-kv11+(kv1pa*F(lm(1))+kv1pb*F(lm(2))+kv1pc*F(lm(3)))+...
(kv1pa*E(lm(1))+kv1pb*E(lm(2))+kv1pc*E(lm(3)))-...
kvpa*(G(lm(1))+kvpb*G(lm(2))+kvpc*G(lm(3)))+...
K*(knn)-M*kv;
k23=sparse(3,3);
k24=sparse(3,3);
k25=lamda*kv;
k26=lamda*N*kv;
%G
k31=sparse(3,3);
k32=sparse(3,3);
k33=sparse(3,3);
k34=sparse(3,3);
k35=sparse(3,3);
k36=sparse(3,3);
%H
k41=sparse(3,3); k42=sparse(3,3);k43=sparse(3,3);
k44=-kv11+kv1pa*F(lm(1))+kv1pb*F(lm(2))+kv1pc*F(lm(3))+...
kv1pa*E(lm(1))+kv1pb*E(lm(2))+kv1pc*E(lm(3))-...
K*(knn)-kvpa*Q(lm(1))+kvpb*Q(lm(2))+kvpc*Q(lm(3))-...
M*kv;
k45=sparse(3,3); k46=sparse(3,3);
%H
k51=sparse(3,3);k52=sparse(3,3);k53=sparse(3,3);
k54=sparse(3,3);
k55=-kv11+(Pr*Nb)*(kv1p1a*H(lm(1))+kv1p1b*H(lm(2))+...
kv1p1c*H(lm(3)))+...
(Pr*Nt)*(kv1p1a*T(lm(1))+kv1p1b*T(lm(2))+...
kv1p1c*T(lm(3)))+...
Pr*kv1pa*F(lm(1))+kv1pb*F(lm(2))+kv1pc*F(lm(3))+...
kv1pa*E(lm(1))+kv1pb*E(lm(2))+kv1pc*E(lm(3))-...
(Pr*delta1)*(kv1pp11*F(lm(1))*G(lm(1))+...
kv1pp12*F(lm(1))*G(lm(2))+...
kv1pp13*F(lm(1))*G(lm(3))+...
kv1pp21*F(lm(2))*G(lm(1))+...
kv1pp22*F(lm(2))*G(lm(2))+...
kv1pp23*F(lm(2))*G(lm(3))+...
kv1pp31*F(lm(3))*G(lm(1))+...
kv1pp32*F(lm(3))*G(lm(2))+...
kv1pp33*F(lm(3))*G(lm(3))+...
kv1pp11*F(lm(1))*Q(lm(1))+...
kv1pp12*F(lm(1))*Q(lm(2))+...
kv1pp13*F(lm(1))*Q(lm(3))+...
kv1pp21*F(lm(2))*Q(lm(1))+...
kv1pp22*F(lm(2))*Q(lm(2))+...
kv1pp23*F(lm(2))*Q(lm(3))+...
kv1pp31*F(lm(3))*Q(lm(1))+...
kv1pp32*F(lm(3))*Q(lm(2))+...
kv1pp33*F(lm(3))*Q(lm(3))+...
kv1pp11*E(lm(1))*G(lm(1))+...
kv1pp12*E(lm(1))*G(lm(2))+...
kv1pp13*E(lm(1))*G(lm(3))+...
kv1pp21*E(lm(2))*G(lm(1))+...
kv1pp22*E(lm(2))*G(lm(2))+...
kv1pp23*E(lm(2))*G(lm(3))+...
kv1pp31*E(lm(3))*G(lm(1))+...
kv1pp32*E(lm(3))*G(lm(2))+...
kv1pp33*E(lm(3))*G(lm(3))+...
kv1pp11*Q(lm(1))*E(lm(1))+...
kv1pp12*Q(lm(1))*E(lm(2))+...
kv1pp13*Q(lm(1))*E(lm(3))+...
kv1pp21*Q(lm(2))*E(lm(1))+...
kv1pp22*Q(lm(2))*E(lm(2))+...
kv1pp23*Q(lm(2))*E(lm(3))+...
kv1pp31*Q(lm(3))*E(lm(1))+...
kv1pp32*Q(lm(3))*E(lm(2))+...
kv1pp33*Q(lm(3))*E(lm(3))+...
kv11pp11*F(lm(1))*F(lm(1))+...
kv11pp12*F(lm(1))*F(lm(2))+...
kv11pp13*F(lm(1))*F(lm(3))+...
kv11pp21*F(lm(2))*F(lm(1))+...
kv11pp22*F(lm(2))*F(lm(2))+...
kv11pp23*F(lm(2))*F(lm(3))+...
kv11pp31*F(lm(3))*F(lm(1))+...
kv11pp32*F(lm(3))*F(lm(2))+...
kv11pp33*F(lm(3))*F(lm(3))+...
2*(kv11pp11*F(lm(1))*E(lm(1))+...
kv11pp12*F(lm(1))*E(lm(2))+...
kv11pp13*F(lm(1))*E(lm(3))+...
kv11pp21*F(lm(2))*E(lm(1))+...
kv11pp22*F(lm(2))*E(lm(2))+...
kv11pp23*F(lm(2))*E(lm(3))+...
kv11pp31*F(lm(3))*E(lm(1))+...
kv11pp32*F(lm(3))*E(lm(2))+...
kv11pp33*F(lm(3))*E(lm(3))))+...
(Pr*delta1)*(kv1pp11*E(lm(1))*E(lm(1))+...
kv11pp12*E(lm(1))*E(lm(2))+...
kv11pp13*E(lm(1))*E(lm(3))+...
kv11pp21*E(lm(2))*E(lm(1))+...
kv11pp22*E(lm(2))*E(lm(2))+...
kv11pp23*E(lm(2))*E(lm(3))+...
kv11pp31*E(lm(3))*E(lm(1))+...
kv11pp32*E(lm(3))*E(lm(2))+...
kv11pp33*E(lm(3))*E(lm(3)));
k56=sparse(3,3);
%Z
k61=sparse(3,3); k62=sparse(3,3);
k63=sparse(3,3); k64=sparse(3,3);
k65=kt*(Nt/Nb);
k66=-kv11+Sc*(kv1pa*F(lm(1))+kv1pb*F(lm(2))+kv1pc*F(lm(3))+...
kv1pa*E(lm(1))+kv1pb*E(lm(2))+kv1pc*E(lm(3))-...
Sc*delta1*(kv1pp11*F(lm(1))*G(lm(1))+...
kv1pp12*F(lm(1))*G(lm(2))+...
kv1pp13*F(lm(1))*G(lm(3))+...
kv1pp21*F(lm(2))*G(lm(1))+...
kv1pp22*F(lm(2))*G(lm(2))+...
kv1pp23*F(lm(2))*G(lm(3))+...
kv1pp31*F(lm(3))*G(lm(1))+...
kv1pp32*F(lm(3))*G(lm(2))+...
kv1pp33*F(lm(3))*G(lm(3))))+...
Sc*delta2*(kv1pp11*F(lm(1))*Q(lm(1))+...
kv1pp12*F(lm(1))*Q(lm(2))+...
kv1pp13*F(lm(1))*Q(lm(3))+...
kv1pp21*F(lm(2))*Q(lm(1))+...
kv1pp22*F(lm(2))*Q(lm(2))+...
kv1pp23*F(lm(2))*Q(lm(3))+...
kv1pp31*F(lm(3))*Q(lm(1))+...
kv1pp32*F(lm(3))*Q(lm(2))+...
kv1pp33*F(lm(3))*Q(lm(3)))+...
Sc*delta2*(kv1pp11*E(lm(1))*G(lm(1))+...
kv1pp12*E(lm(1))*G(lm(2))+...
kv1pp13*E(lm(1))*G(lm(3))+...
kv1pp21*E(lm(2))*G(lm(1))+...
kv1pp22*E(lm(2))*G(lm(2))+...
kv1pp23*E(lm(2))*G(lm(3))+...
kv1pp31*E(lm(3))*G(lm(1))+...
kv1pp32*E(lm(3))*G(lm(2))+...
kv1pp33*E(lm(3))*G(lm(3)))+...
Sc*delta2*(kv1pp11*Q(lm(1))*E(lm(1))+...
kv1pp12*Q(lm(1))*E(lm(2))+...
kv1pp13*Q(lm(1))*E(lm(3))+...
kv1pp21*Q(lm(2))*E(lm(1))+...
kv1pp22*Q(lm(2))*E(lm(2))+...
kv1pp23*Q(lm(2))*E(lm(3))+...
kv1pp31*Q(lm(3))*E(lm(1))+...
kv1pp32*Q(lm(3))*E(lm(2))+...
kv1pp33*Q(lm(3))*E(lm(3)))+...
Sc*delta2*(kv11pp11*F(lm(1))*F(lm(1))+...
kv11pp12*F(lm(1))*F(lm(2))+...
kv11pp13*F(lm(1))*F(lm(3))+...
kv11pp21*F(lm(2))*F(lm(1))+...
kv11pp22*F(lm(2))*F(lm(2))+...
kv11pp23*F(lm(2))*F(lm(3))+...
kv11pp31*F(lm(3))*F(lm(1))+...
kv11pp32*F(lm(3))*F(lm(2))+...
kv11pp33*F(lm(3))*F(lm(3)))+...
2*(kv11pp11*F(lm(1))*E(lm(1))+...
kv11pp12*F(lm(1))*E(lm(2))+...
kv11pp13*F(lm(1))*E(lm(3))+...
kv11pp21*F(lm(2))*E(lm(1))+...
kv11pp22*F(lm(2))*E(lm(2))+...
kv11pp23*F(lm(2))*E(lm(3))+...
kv11pp31*F(lm(3))*E(lm(1))+...
kv11pp32*F(lm(3))*E(lm(2))+...
kv11pp33*F(lm(3))*E(lm(3)))+...
Sc*delta2*(kv1pp11*E(lm(1))*E(lm(1))+...
kv11pp12*E(lm(1))*E(lm(2))+...
kv11pp13*E(lm(1))*E(lm(3))+...
kv11pp21*E(lm(2))*E(lm(1))+...
kv11pp22*E(lm(2))*E(lm(2))+...
kv11pp23*G(lm(2))*E(lm(3))+...
kv11pp31*E(lm(3))*E(lm(1))+...
kv11pp32*E(lm(3))*E(lm(2))+...
kv11pp33*E(lm(3))*E(lm(3)));
l1(lm,lm)=l1(lm,lm)+k11;
l2(lm,lm)=l2(lm,lm)+k12;
l3(lm,lm)=l3(lm,lm)+k13;
l4(lm,lm)=l4(lm,lm)+k14;
l5(lm,lm)=l5(lm,lm)+k15;
l6(lm,lm)=l6(lm,lm)+k16;
l7(lm,lm)=l7(lm,lm)+k21;
l8(lm,lm)=l8(lm,lm)+k22;
l9(lm,lm)=l9(lm,lm)+k23;
l10(lm,lm)=l10(lm,lm)+k24;
l11(lm,lm)=l11(lm,lm)+k25;
l12(lm,lm)=l12(lm,lm)+k26;
l13(lm,lm)=l13(lm,lm)+k31;
l14(lm,lm)=l14(lm,lm)+k32;
l15(lm,lm)=l15(lm,lm)+k33;
l16(lm,lm)=l16(lm,lm)+k34;
l17(lm,lm)=l17(lm,lm)+k35;
l18(lm,lm)=l18(lm,lm)+k36;
l16(lm,lm)=l19(lm,lm)+k41;
l20(lm,lm)=l20(lm,lm)+k42;
l21(lm,lm)=l21(lm,lm)+k43;
l22(lm,lm)=l22(lm,lm)+k44;
l23(lm,lm)=l23(lm,lm)+k45;
l24(lm,lm)=l24(lm,lm)+k46;
l25(lm,lm)=l25(lm,lm)+k51;
l26(lm,lm)=l26(lm,lm)+k52;
l27(lm,lm)=l27(lm,lm)+k53;
l28(lm,lm)=l28(lm,lm)+k54;
l29(lm,lm)=l29(lm,lm)+k55;
l30(lm,lm)=l30(lm,lm)+k56;
l31(lm,lm)=l31(lm,lm)+k61;
l32(lm,lm)=l32(lm,lm)+k62;
l33(lm,lm)=l33(lm,lm)+k63;
l34(lm,lm)=l34(lm,lm)+k64;
l35(lm,lm)=l35(lm,lm)+k65;
l36(lm,lm)=l36(lm,lm)+k66;
R1(lm)=R1(lm)+ f1;
R2(lm)=R2(lm)+ f2;
R3(lm)=R3(lm)+ f3;
R4(lm)=R4(lm)+ f4;
R5(lm)=R5(lm)+ f5;
R6(lm)=R6(lm)+ f6;
end
l=[l1,l2,l3,l4,l5,l6;l7,l8,l9,l10,l11,l12;l13,l14,l15,l16,l17,l18;l19,l20,l21,l22,l23,l24;l25,l26,l27,l28,l29,l30;l31,l32,l33,l34,l35,l36];
LH=zeros(14,6*nn);
RH=zeros(6*nn,14);
RC=zeros(14,14);
LH(1,1)=1;
LH(2,nn+1)=1;
LH(3,2*nn)=1;
LH(4,2*nn+1)=1;
LH(5,3*nn)=0;
LH(6,3*nn+1)=1;
LH(7,4*nn)=1;
LH(8,4*nn+1)=1;
LH(9,5*nn)=1;
LH(10,5*nn+1)=1;
LH(11,6*nn)=1;
LH(12,6*nn+1)=1;
LH(13,7*nn)=1;
LH(14,7*nn+1)=1;
RH(1,1)=1;
RH(nn+1,2)=-1;
RH(2*nn,3)=1;
RH(2*nn+1,4)=-1;
RH(3*nn,5)=1;
RH(3*nn+1,6)=-1;
RH(4*nn,7)=1;
RH(4*nn+1,8)=-1;
RH(5*nn,9)=1;
RH(5*nn+1,10)=-1;
RH(6*nn,11)=1;
RH(6*nn+1,12)=-1;
RH(7*nn,13)=1;
RH(7*nn+1,14)=-1;
%R=[R1;R2;R3;R4;R5;R6;0;1;0;(beta1*G(nn+1));0;(1+(T(nn+1)));0;(1+(H(nn+1)));0;0;0];
%R=[R1;R2;R3;R4;R5;R6;0;1;0;1;0;0;-K*G(nn+1);K*G(nn+1);1;0;0];Nb*H(nn+1)+Nt*T(nn+1)
R=[R1;R2;R3;R4;R5;R6;0;1+gamma*G(nn+1);0;0;0;0;(alpha+beta1*Q(nn+1));0;0;0;(1+(T(nn+1)/Bi));0;-(Nt/Nb)*T(nn+1);0];
% R=[R1;R2;R3;R4;R5;R6;0;1+gamma*G(nn+1);0;(alpha+beta1*Q(nn+1));0;(1+(T(nn+1)/Bi));0;1;0;1;(1+(H(nn+1)/Bi))];
%R=[R1;R2;R3;R4;R5;R6;0;(1+gamma*G(nn+1));0;(alpha+beta1*Q(nn+1));0;1;0;1;0;0;0];
l=[l,RH;LH,RC];
L=sparse(l);
% SOLVE FOR NODAL PARAMETERS
d = L\R;
F1 = [d(1:nn);d(6*nn+1);d(6*nn+2)];
G1 = [d(nn+1:2*nn);d(6*nn+3);d(6*nn+4)];
E1 = [d(2*nn+1:3*nn);d(6*nn+5);d(6*nn+6)];
Q1 = [d(3*nn+1:4*nn);d(6*nn+7);d(6*nn+8)];
T1 = [d(4*nn+1:5*nn);d(6*nn+19);d(6*nn+10);d(6*nn+11)];
H1 = [d(5*nn+1:6*nn);d(6*nn+12);d(6*nn+13);d(6*nn+14)];
end
function [kv11,kv1,kv,kv1pa,kv1pb,kv1pc,kvpa,kvpb,kvpc,...
kv1pp11,kv1pp12,kv1pp13,kv1pp21,kv1pp22,kv1pp23,kv1pp31,...
kv1pp32,kv1pp33,kv11pp11,kv11pp12,kv11pp13,kv11pp21,kv11pp22,...
kv11pp23,kv11pp31,kv11pp32,kv11pp33,...
kv1p1a,kv1p1b,kv1p1c,kt,knn] = quadraticelement(he,i)
gL = [-sqrt(3/5);0;sqrt(3/5)];
gW = [5/9,8/9,5/9 ];
kv11=sparse(3,3);
kv1=sparse(3,3);
kv=sparse(3,3);
kv1pa=sparse(3,3);
kv1pb=sparse(3,3);
kv1pc=sparse(3,3);
kvpa=sparse(3,3);
kvpb=sparse(3,3);
kvpc=sparse(3,3);
kv1pp11=sparse(3,3);
kv1pp12=sparse(3,3);
kv1pp13=sparse(3,3);
kv1pp21=sparse(3,3);
kv1pp22=sparse(3,3);
kv1pp23=sparse(3,3);
kv1pp31=sparse(3,3);
kv1pp32=sparse(3,3);
kv1pp33=sparse(3,3);
kv11pp11=sparse(3,3);
kv11pp12=sparse(3,3);
kv11pp13=sparse(3,3);
kv11pp21=sparse(3,3);
kv11pp22=sparse(3,3);
kv11pp23=sparse(3,3);
kv11pp31=sparse(3,3);
kv11pp32=sparse(3,3);
kv11pp33=sparse(3,3);
kv1p1a=sparse(3,3);
kv1p1b=sparse(3,3);
kv1p1c=sparse(3,3);
kt=sparse(3,3);
knn=sparse(3,3);
for k=1:length(gW)
s = gL(k); w = gW(k);
n = [-(1 - s)*s/2,(1-s^2),(1 + s)*s/2];
dns=[(2*s-1)/2 ,-2*s,(2*s+1)/2];
ddns=[1,-2,1];
coord=[(i-1)*he;(i-1/2)*he;i*he];
x = n*coord(:,1);
J=he/2;
dx=dns*(1/J);
ddx=ddns*(1/J)*(1/J);
kv11=kv11+J*w*dx'*dx;
kv1=kv1+J*w*n'*dx;
kv=kv+J*w*n'*n;
kv1pa=kv1pa+J*w*n'*dx*n(1);
kv1pb=kv1pb+J*w*n'*dx*n(2);
kv1pc=kv1pc+J*w*n'*dx*n(3);
kvpa=kvpa+J*w*n'*n*n(1);
kvpb=kvpb+J*w*n'*n*n(2);
kvpc=kvpc+J*w*n'*n*n(3);
kv1pp11=kv1pp11+J*w*n'*dx*n(1)*n(1);
kv1pp12=kv1pp12+J*w*n'*dx*n(1)*n(2);
kv1pp13=kv1pp13+J*w*n'*dx*n(1)*n(3);
kv1pp21=kv1pp21+J*w*n'*dx*n(2)*n(1);
kv1pp22=kv1pp22+J*w*n'*dx*n(2)*n(2);
kv1pp23=kv1pp23+J*w*n'*dx*n(2)*n(3);
kv1pp31=kv1pp31+J*w*n'*dx*n(3)*n(1);
kv1pp32=kv1pp32+J*w*n'*dx*n(3)*n(2);
kv1pp33=kv1pp33+J*w*n'*dx*n(3)*n(3);
kv11pp11=kv11pp11+J*w*n'*ddx*n(1)*n(1);
kv11pp12=kv11pp12+J*w*n'*ddx*n(1)*n(2);
kv11pp13=kv11pp13+J*w*n'*ddx*n(1)*n(3);
kv11pp21=kv11pp21+J*w*n'*ddx*n(2)*n(1);
kv11pp22=kv11pp22+J*w*n'*ddx*n(2)*n(2);
kv11pp23=kv11pp23+J*w*n'*ddx*n(2)*n(3);
kv11pp31=kv11pp31+J*w*n'*ddx*n(3)*n(1);
kv11pp32=kv11pp32+J*w*n'*ddx*n(3)*n(2);
kv11pp33=kv11pp33+J*w*n'*ddx*n(3)*n(3);
kv1p1a=kv1p1a+J*w*n'*dx*dx(1);
kv1p1b=kv1p1b+J*w*n'*dx*dx(2);
kv1p1c=kv1p1c+J*w*n'*dx*dx(3);
kt=kt+J*w*n'*ddx;
knn=knn+J*w*ddx'*ddx;
end
end
%--End
Where is my error when I run math give this
Dimensions of arrays being concatenated are
not consistent.
Error in couple130m>assembly (line 488)
l=[l,RH;LH,RC];
Error in couple130m (line 42)
[F1,G1,E1,Q1,T1,H1] =
assembly(F,G,E,Q,T,H,n,nn,he,Pr,M,Sc,gamma,...
Torsten
Torsten el 5 de Nov. de 2023
Movida: John D'Errico el 5 de Nov. de 2023
At the moment, there can be two reasons if your code doesn't work: your mathematical problem is hard to solve or ill-posed or your finite-element code has errors.
So I repeat my advice:
Start setting up your problem first within "bvp4c" or "bvp5c". If you get problems or errors here, you can be quite sure that they are due to your problem formulation and not due to the method used to solve it. Then - at the time you have a working model within "bvp4c" or "bvp5c" - you might dare to program your own solver.

Iniciar sesión para comentar.

Categorías

Más información sobre Boundary Value Problems en Help Center y File Exchange.

Etiquetas

Productos

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by