what is finite elemnt method code for Sytem of ODE

2 views (last 30 days)
Tadese
Tadese on 22 Oct 2023
Moved: John D'Errico on 5 Nov 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.

Sign in to answer this question.

Answers (1)

Torsten
Torsten on 22 Oct 2023
Moved: Torsten on 22 Oct 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 Comments
Show 2 older comments
Tadese
Tadese on 5 Nov 2023
Moved: John D'Errico on 5 Nov 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 on 5 Nov 2023
Moved: John D'Errico on 5 Nov 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.

Sign in to comment.

Categories

Find more on Boundary Value Problems in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!