what is finite elemnt method code for Sytem of ODE

2 次查看（过去 30 天）
Tadese
Tadese 2023-10-22
移动：John D'Errico 2023-11-5
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.

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2023-10-22
移动：Torsten 2023-10-22
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 个评论
显示 2更早的评论
Tadese
Tadese 2023-11-5
移动：John D'Errico 2023-11-5
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 2023-11-5
移动：John D'Errico 2023-11-5
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by