Comparison PDE solve and \ mldivide

1 次查看（过去 30 天）
Phineas
Phineas 2025-9-29
编辑： Walter Roberson 2025-9-29
Ran in:
Hello all,
I have a problem and would like you thoughts on it. I guess I miss something in the definition of some matlab functions but can not find out where my error is.
Comparison of solve function with mldivide in a transient heat analysis gives different solutions.
Probably an missuse of the Nullspace transformation but did not find yet where??
I would really like to understand what I got wrong and would be thankfull for your comments.
Phineas Freak
%Following unified modelling approach:
clear all
%Geometry
shape_1=polyshape([-150, 150, 150, -150],[-150, -150, 150,150]);
BlockRaw= triangulation(shape_1);
BlockRaw=fegeometry(BlockRaw);
BlockRaw=extrude(BlockRaw,300);
BlockRaw=scale(BlockRaw, 1/1000); % [mm] to [m]
% Model
model=femodel(AnalysisType="ThermalTransient", Geometry= BlockRaw);
figure
pdegplot(model,FaceLabels="on",EdgeLabels="on")
model.MaterialProperties=materialProperties(ThermalConductivity=46.5, MassDensity=7800, SpecificHeat=.502);
% Initial Temperature
model.CellIC(1) = cellIC(Temperature=273+180);
%BC1
model.FaceLoad(3) = faceLoad(Heat=300);
% BC2
model.FaceBC(1) = faceBC(Temperature=200);
% BC 3
model.FaceLoad(6) = faceLoad(Heat=100);
% Coarse Mesh - exact result not important as comparison between methods
model = generateMesh(model, GeometricOrder='linear', Hmin=.05);
figure
pdemesh(model)
% Time Interval
dt=.25;
t_end=4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Builtin solver
ZylinderRes = solve(model, [0:dt:t_end]);
figure
pdeplot3D(ZylinderRes.Mesh,ColorMapData=ZylinderRes.Temperature(:,end))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mldivide
state = struct('time', ZylinderRes.SolutionTimes(1), 'u', ZylinderRes.Temperature(:,1)');
FEMMatrix=assembleFEMatrices(model,state);
K = FEMMatrix.K;
M= FEMMatrix.M;
F = FEMMatrix.F;
H=FEMMatrix.H;
R=FEMMatrix.R;
G=FEMMatrix.G;
theta=1;
[B,Or]=pdenullorth(FEMMatrix.H);
ud=Or*((H*Or\R)); % Prüfung - Korrekt
U_=zeros(length(ud),17);
U_(:,1)=ZylinderRes.Temperature(:,1);
% Nullspace Dirichlet Transformation
Kc =B'*(theta*K + M/dt)*B;
dt_sum = 0;
for i=2:t_end/dt+1
dt_sum=dt_sum+dt;
Fc=B'*((F+G) - K*ud + ( M/dt - (1-theta)*K )*U_(:,i-1)) ;
U_(:,i) = B* (Kc \ Fc ) + ud;
end
figure
pdeplot3D(ZylinderRes.Mesh,ColorMapData=U_(:,end))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Comparison: Missmatch
figure
plot(ZylinderRes.Temperature(1:400,4)-U_(1:400,4),'b-')
  3 个评论
显示 1更早的评论
Phineas
Phineas 2025-9-29
编辑：Phineas 2025-9-29
Hello Torsten,
Thanks for the redaction of the question, looks smooth.
  1. I can not edit the solve function, its a MathWorks built in function. I took the proceeding for solve from the examples.
  2. Calculation of U_ follows different literature sources that are quite coherent (except me making an error). [1] , [2] , [3].
Problem here is, I do not know how M from assembleFEMatrices() is calculated, Matlab manual is somehow ambiguous here. Can be coefficient for d2T/dt2 or dT/dt term: "M is the mass matrix, the integral of the discretized version of the m or d coefficients. " --> Section ALGORITHMS
Phineas
Phineas 2025-9-29
In fact, I will try solvepde() might be this changes the solution....

请先登录，再进行评论。

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

采纳的回答

Phineas
Phineas 2025-9-29
编辑：Walter Roberson 2025-9-29
This seems to work.
clear all
dt=.25;
t_end=4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
model2=createpde();
model2.Geometry=multicuboid(.3,.3,.3);
pdegplot(model2.Geometry,FaceLabels="on",FaceAlpha=0.5)
axis equal
model2.generateMesh(GeometricOrder='linear',Hmin=.05);
setInitialConditions(model2,273+180);
applyBoundaryCondition(model2,"dirichlet", ...
"Face",1,"u",100)
applyBoundaryCondition(model2,"neumann",Face=[3],q=0,g=200);
specifyCoefficients(model2,m=0,...
d=1,...
c=1,...
a=0,...
f=1);
Result2 = solvepde(model2, [0:dt:t_end]);
pdeplot3D(model2,"ColorMapData",Result2.NodalSolution(:,5))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Own Solver
FEMMatrix2=assembleFEMatrices(model2)
K = FEMMatrix2.K;
M = FEMMatrix2.M;
H = FEMMatrix2.H;
R = FEMMatrix2.R;
G = FEMMatrix2.G;
%theta=1;
[B,Or]=pdenullorth(H);
ud=Or*((H*Or\R));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
U_=zeros(length(ud),t_end/dt+1);
U_(:,1)=Result2.NodalSolution(:,1);
theta=1;
Kc =B'*(theta*K + M/dt)*B;
dt_sum=0;
for i=2:t_end/dt+1
dt_sum=dt_sum+dt;
Fc=B'*( G + ( M/dt -(1-theta)*K )*U_(:,i-1) - K*ud) ;
U_(:,i) = B* (Kc \ Fc ) + ud;
end
for i=1:t_end/dt+1
mean(Result2.NodalSolution(:,i)-U_(:,i))
end
for i=1:t_end/dt+1
figure
plot(Result2.NodalSolution(1:400,10),U_(1:400,10),'b*')
end

更多回答（1 个）

Walter Roberson
Walter Roberson 2025-9-29
model.MaterialProperties(ThermalConductivity=46.5, MassDensity=7800, SpecificHeat=.502);
That is invalid. It is converted internally into
model.MaterialProperties("ThermalConductivity", 46.5, "MassDensity", 7800, "SpecificHeat", .502);
which attempts to use the string "ThermalConductivity" as an index into the model.MaterialProperties array.
You instead need something like
model.MaterialProperties = materialProperties(ThermalConductivity=46.5, ...
MassDensity=7800, ...
SpecificHeat=0.502);
  3 个评论
显示 1更早的评论
Phineas
Phineas 2025-9-29
I changed it but for the moment still the gap. Will switch doe solvepde().
Walter Roberson
Walter Roberson 2025-9-29
The call you had was simply invalid. The replacement code uses the properties that you were trying to specify, not different properties.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Heat Transfer 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by