- Are you looking into the equation1) rho * Cp * (du/dt) = k * (d^2u/dt^2) + s? and 2) deltaT/deltat = (s*L) / (h * A)?
- If we are rewriting the equation 2) like(Potf - Poti) / (tf - ti) = (s * L) / (h * A)then with the available value of Potf, Poti, tf, ti, h and A, you can find the value of s right?
You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
I need to know how to put the s variable of pdefun (from pdepe solver) as a function of x and time
2 views (last 30 days)
Show older comments
Hi, I need help to solve a proble I'm in for a long ago. The thing is that I can't find any doc explaining how to set the variables [c,f,s] of pdefun (in particular the 's') as a function of the point of the 'xmesh'. It is to solve the heat transfer and difussion equation that is rho*cp*du/dt-k*d/dx(du/dx)=g.
from mathworks blogs and answers, I know that c=rho*cp (which is constant with time and x) f=k (also constant), but g is power generation, and that is not constant with time neither with x, the power is generated the first time values and just in a given set of points of xmesh.
Here I paste a part of the code where I state the pdefun:
function[c,f,s] = pdef(x, u, t, DuDx)
c = rho*cp;
f = k*DuDx;
i=1;
while(i<=length(x))
if(x(i)<=percent*radi)
s = (1/(h*pi*radi)^2)).*interp1(t40,Pot40,t,'linear','extrap');
else
s = 0;
end
i=i+1;
end
end
Please, if anyone could tell me if this is well developed. The s variable I mean.
8 Comments
Nikhil Sreekumar
on 9 Jan 2017
Hi Joan,
I had some doubts could you please clarify?
Thanks
Nikhil
JOAN PERE PONSETI
on 10 Jan 2017
Hi, yes, my doubt is regarding to the equation you call 1), but variable u in my case is temperature T. The thing is that, matlab says that s is a function of x,u,t,dudx and I just want it to put s(x,t) because the source term depends of the position within the radius (x) and the time t, with the power. My source term is the power density generation, which is P/volume and the power depends on time and the volume depends on x. this is why my source term depends on x and time.
thanks!
JP
Torsten
on 11 Jan 2017
The formulation
"s is a function of x,u,t,dudx"
means that s may depend on all of these variables, but of course it is possible that it depends on x and t only.
Best wishes
Torsten.
JOAN PERE PONSETI
on 11 Jan 2017
Yes, I know, but my problem is that I don't know exactly how the anonymous function of pdefun works, I mean, if the pdepe calls for the pdefun each time in 'tspan' and each x in 'xmesh', if so, pdefun returns 's' as a single value (double) but don't know how to give the function the vectors of tspan and xmesh and then return that single value.
Thanks!
JOAN PERE PONSETI
on 11 Jan 2017
The thing is that there is a velue for s in each time and each x, then, for a given x point there will be never the same value of s because it also changes on time. So, if I have a tspan of 20 intervals, and and xmesh of 50 points, s should be a matrix (50x20) which will have in each cell the value of the source term at each time and each x.
Torsten
on 11 Jan 2017
Edited: Torsten
on 11 Jan 2017
Why do you think you have to supply tspan and xmesh in "pdefun" ?
From the documentation of pdepe:
pdefun computes the terms c, f, and s (Equation 1-3). It has the form
[c,f,s] = pdefun(x,t,u,dudx)
The input arguments are scalars x and t and vectors u and dudx that approximate the solution u and its partial derivative with respect to x, respectively. c, f, and s are column vectors. c stores the diagonal elements of the matrix c (Equation 1-3).
So your code should read:
function[c,f,s] = pdef(x, u, t, DuDx)
c = rho*cp;
f = k*DuDx;
if x<=percent*radi
s = (1/(h*pi*radi)^2))*interp1(t40,Pot40,t,'linear','extrap');
else
s = 0;
end
end
Of course, the variables percent, radi, h, t40 and Pot40 must be visible in "pdef".
Best wishes
Torsten.
JOAN PERE PONSETI
on 11 Jan 2017
Hi again, I implemented this code to see if now it works but it gives me this error message back.
Not enough input arguments.
Error in CLDissipation1D/pdef (line 375) if(x==xmesh(1))
Error in pdepe (line 246) [c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});
Error in CLDissipation1D (line 340) sol = pdepe(m, @pdef, @pdeic, @pdebc, xmesh, t40);
function[c,f,s] = pdef(x,t,DuDx,S,xmesh,tspan)
c=rho*cp;
f=k*DuDx;
% find out the volume we need the previous point too
i=1;
encontrado=0;
while(encontrado==0)
if(x==xmesh(1))
vol=inf;
encontrado=1;
elseif(x==xmesh(i))
vol=(pi*(x^2-xmesh(i-1)^2)*h);
encontrado=1;
end
i=i+1;
end
% repeat the same loop to find the power at time t
i=1;
encontrado=0;
while(encontrado==0)
if(t==tspan(i))
pot=Pot40(i);
end
i=i+1;
end
% Finally we define the variable s as the ratio of power/volum
s=pot/(vol*pin);
% ========================================================================= % defines the Initial Conditions
end
Thanks a lot again!
Accepted Answer
Torsten
on 11 Jan 2017
1. pdepe calls pdef for arbitrary times, not only those specified in tspan.
2. sol = pdepe(m, @(x,t,DuDx)pdef(x,t,DuDx,S,xmesh,tspan), @pdeic, @pdebc);
3. Comparisons for equality like if x==xmesh(1) probably won't work because of rounding errors.
Best wishes
Torsten.
9 Comments
JOAN PERE PONSETI
on 13 Jan 2017
Ok, I have tried this, but matlab returns this error.
Not enough input arguments.
Error in CLDissipation1D>@(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh) (line 296) sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);
Error in pdepe (line 246) [c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});
Error in CLDissipation1D (line 296) sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);
sol = pdepe(m, @(x,t,DuDx,vol,xmesh)pdef(x,t,DuDx,vol,xmesh), @pdeic, @pdebc, xmesh, t40);
function[c,f,s] = pdef(x,t,DuDx,vol,xmesh)
c=rho*cp;
f=k*DuDx;
s=interp1(t40,Pot40,t,'linear','extrap')/(pin.*interp1(xmesh,vol,x,'linear','extrap'));
% ========================================================================= % defines the Initial Conditions
end
Thanks a lot!
JOAN PERE PONSETI
on 13 Jan 2017
And now, this. Here the problem is in pdeic??
Error using CLDissipation1D/pdeic Too many input arguments.
Error in pdepe (line 229) temp = feval(ic,xmesh(1),varargin{:});
Error in CLDissipation1D (line 188) sol = pdepe(m, @(vol)pdef,@pdeic, @pdebc, xmesh, t40,[], vol);
sol = pdepe(m, @(vol)pdef,@pdeic, @pdebc, xmesh, t40,[], vol);
u=sol(:,:,1);
%========================================================================= % defines the properties of the PDE function
function[c,f,s] = pdef(x,t,u,DuDx,vol)
c=rho*cp;
f=k*DuDx;
s=interp1(t40,Pot40,t,'linear','extrap')/(pin.*interp1(xmesh,vol,x,'spline','extrap'));
% ========================================================================= % defines the Initial Conditions
end
function[u0] = pdeic(x)
u0 = tinitial;
end % ========================================================================= % defines the Boundary Conditions
function[pl,ql,pr,qr] = pdebc(xl, ul, xr, ur, t)
%for cylindrical pl and ql are 0, because of symmetry. No flux is present
%at the center of the sphere.
pl = 0;
ql = 0;
pr = ur-interp1(t40,Temp40,t,'linear','extrap');
qr = 0;
end
Torsten
on 16 Jan 2017
Take a look at this page on how to pass extra parameters by the anonymous function method:
https://de.mathworks.com/help/optim/ug/passing-extra-parameters.html
Best wishes
Torsten.
JOAN PERE PONSETI
on 21 Jan 2017
Thanks a lot! I've solved this by stating the variables as global. But now, I have another problem with the solution I get. The thing is that, there is a power dissipation at the center, and I have the initial conditions and the boundary ones. Logic says that if there is power dissipation inside, the temperature within the core must be always higher than the one outside (at the same time for both temperatures). But what I get is the graph below, and furthermore, the blue-crossed line of the graph doesn't start at the 300 K (aprox) which is the initial condition.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/185519/image.png)
I give you the code I have used and let's see if you can figure out what's the error.
Thanks!
Joan Pere
It looks like the differential is not working properly.. don't know ehat happens
t40=times{var};
Temp40=temperatures{var};
Pot40=potencies{var};
% Assignem els vectors a variables que es mostraran en el Workspace
assignin('base','t40',t40);
assignin('base','Temp40',Temp40);
assignin('base','Pot40',Pot40)
% figure
% plot(t40,Pot40);
%
% figure
% plot(t40,Temp40);
% capacitor properties
k = 1.5; % W/(m K) 30 Alumina %aluminum 237
rho = 3000; % alumina density 4000 Aluminum 2300
cp = 900; % J/K/kg Aluminum 910 Alumina 880
%capacitor geometry
h=0.01;
radi=0.005;
percent=0.25;
tinitial = Temp40(1); % intial temp throughout the tile
m = 1; % assume cylindrical symmetry
xmesh = linspace(0, radi, 200); % 100 punts a discretitzar
assignin('base','xmesh',xmesh)
global vol
for i=1:length(xmesh)
if(i==1)
vol(i)=pi*h*(radi)^2;
else
vol(i)=inf;
end
end
assignin('base','vol',vol)
sol = pdepe(m, @pdef,@pdeic, @pdebc, xmesh, t40);
T1=sol(:,:,1);
assignin('base','T1',T1)
% A surface plot is often a good way to study a solution.
figure;
subplot(2,1,1)
surf(xmesh,t40,T1);
title(sprintf('SOLUCIÓ NUMÈRICA AMB 200 PUNTS.\nArxiu %s.', VCC40));
xlabel('Distància al centre (m)');
ylabel('Temps (segons)');
% A solution profile can also be illuminating.
subplot(2,1,2)
plot(t40,T1(:,1),'x');
hold on
plot(t40,T1(:,end),'+');
hold on
plot(t40,Temp40,'o');
legend('Temperatura al centre','Temperatura externa','Condicions Incials','Location','best');
title(sprintf('PERFIL DE TEMPERATURES EN 2 PUNTS (x=0 i x=R).\nArxiu %s.', VCC40));
xlabel('Temps (segons)');
ylabel('Temperatura (Kelvin)');
%========================================================================== % defines the properties of the PDE function
function[c,f,s] = pdef(x,u,t,DuDx)
c=rho*cp;
f=k*DuDx;
s=(1/interp1(xmesh,vol,x,'linear','extrap')).*interp1(t40,Pot40,t,'linear','extrap');
end
% ========================================================================= % defines the Initial Conditions
function[u0] = pdeic(x)
u0 = tinitial;
end
% ========================================================================= % defines the Boundary Conditions
function[pl,ql,pr,qr] = pdebc(xl, ul, xr, ur, t)
%for cylindrical pl and ql are 0, because of symmetry. No flux is present at the center of the sphere.
pl = 0;
ql = 0;
pr = ur-interp1(t40,Temp40,t,'linear','extrap');
qr = 0;
end
The volume is needed because my source term is the power density (Power/Volume) then the power is a different one each time (that's why I have to interpolate) and the volume I state that the volume dissipation is a quarter of the radius (percent=0.25)
Torsten
on 23 Jan 2017
1. In the boundary part (@pdebc), pl=ql=0 can't be correct.
2. I don't see where "percent" is used in the calculation of the source term.
3. I don't understand your variable "vol" because you assign "inf" to all the inner grid points.
Best wishes
Torsten.
JOAN PERE PONSETI
on 24 Feb 2017
1. pl and ql must be zero because there is no te flux into that point, which is the center of the cilindre.
2. the percent variable let's say is irrelevant here.
3. I see now that the volume variable code is not complete, I might have deleted some code lines. But the thing is that, i need the inner points to have a given value (power/volume, both of them are scalars) but i need the outter ones to be 0, thats why I set de volume to inf, since it divides and makes the value of the source term 0 (power always have a positive value).
Thanks!
JOAN PERE PONSETI
on 24 Feb 2017
Edited: JOAN PERE PONSETI
on 24 Feb 2017
Hello again, I need some help in another code, similar to this one but now using the solvepde function. When I have to assign the coefficients. I guess that when I write the following code:
sourcefun= @(region,state) interp1(t40,Pot40,state.t,'linear','extrap')/(volumeD*count);
specifyCoefficients(pdem,'m',0,'d',rho*cp,'c',k,'a',0,'f', sourcefun);
Matlab returs an error that says the following:
Error using pde.CoefficientAssignment/checkCoefFcnHdlArgCounts (line 531)
Function handle specifying a coefficient must accept two
input arguments and return one output argument.
As if the function doesn't return any value, is it parmited to use the state and region variables in a function like interp1??
Thanks a lot!!
Joan Pere
JOAN PERE PONSETI
on 24 Feb 2017
Another thing, when I write the next code % Specify Coefficients
f = @(region,state) interp1(t40,Pot40,state.time,'linear','extrap');
specifyCoefficients(pdem,'m',0,'d',rho*cp,'c',k,'a',0,'f', f);
and I set a breakpoint at f, it always says that state structure is 0 at all its elements but it doesn't return any error until the fifth or sixth step... don't know why...
thanks!
JClarcq
on 19 Feb 2018
Hi Joan,
I have the same problem than you where my heatsource is a term temperature dependent via an interp1 of a material properties.
Have you been able to solve this inputs argument issue? Thanks,
More Answers (0)
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)