Laplace Equation with pdepe command

Hi, I am learning to solve PDEs with pdepe command. First, I try to solve a Laplace equation (1D) setting zero the right coefficients in the PDE. I am expecting the solution to be e.g. y=m*x + b. However I get curves and not lines and I do not know why. I attach the M-files below.

Thanks.

Answer 1

[Torsten]

Setting s=1, the solution is of the form u(x)=-1/2*x^2+a*x+b.

Best wishes

Torsten.

[Thanasis Hou]

I corrected it but I did not get the desire result. I made a simulation with Comsol(I attached it). Could I get the same result with the MatLab using the pdepe command, making a Laplace Equation.

[Torsten]

It's written in the documentation that c identically zero is not allowed. So you could try adding a second (artificial) equation (e.g. du2/dt = 0).

From the documentation:

In Equation 1-3, f (x,t,u,∂u/∂x) is a flux term and s (x,t,u,∂u/∂x) is a source term. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix c (x,t,u,∂u/∂x). The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if those values of x are mesh points. Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface.

Best wishes

Torsten.

[Torsten]

If your equation is simply elliptic, use "bvp4c" instead of "pdepe".

Best wishes

Torsten.

[Thanasis Hou]

Thank you very much. I will look better on the Internet.

Thanasis

[MEENAL SINGHAL]

@Torsten

Your suggestion "It's written in the documentation that c identically zero is not allowed. So you could try adding a second (artificial) equation (e.g. du2/dt = 0)." seems valuable to me.

I was wondering what would be the initial condition, if i introduce du2/dt = 0, in addition to my previous equations such that the overall system (a steady state heat transfer equation) remains the same?

Thank you!

-Meenal

[MEENAL SINGHAL]

To add on the complexity, my equations contain discontinuity.

Pdes are defined as in pdes.jpg and the corresponding code is

function forward_compositewalls

global Nci Nco n Nri Nro alphai Betai Gammai alphao Betao Gammao Deltai Deltao et_i et_o 
global Xbw Qi Qo theta_c theta_rad theta_ini Temperature_1 Temperature

Xbw=0.6;
Nci=0; Nco=1;
n=0.5;
Nri=1; Nro=1;
alphai=0.01; alphao=0.01;
Betai=0.04; Betao=0.04;
Gammai=0.05; Gammao=0.03;
Deltai=0.01; Deltao=0.01;
%kratio=0.01;
Qi=1; Qo=1;
theta_c=0.025;
theta_rad=0.025; 
theta_ini=0.03;
et_i=0.6;
et_o=0.4;
L=1;
tend = 3;
m=0;

x = linspace(0,L,21);

t= linspace(0,tend,11);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);


% Extract the solution components as Temperature.
Temperature = sol(:,:,1);
Temperature_1=Temperature';
% A solution profile can also be illuminating.
figure, plot(x,(Temperature_1(:,end)))
% 

function [c,f,s] = pdex1pde(x,t,u,DuDx)
global alphai Betai Gammai alphao Betao Gammao Deltai Deltao Xbw Qi Qo theta_rad 

if x<=Xbw
    c = [0;1];
    f = [1 + Deltai*(u-1).*(DuDx);0];
    s = [Qi.*(1+alphai*u+Betai*u.^2+Gammai*u.^3);0];
else
    c = [0;1];
    f = [1 + Deltao*(u-theta_rad).*(DuDx);0];
    s = [Qo*(1+alphao*u+Betao*u.^2+Gammao*u.^3);0];
end
function u0 = pdex1ic(x)

u0 =[ 0; 0 ];  % i don't know this (What should be the value for the steady case?)
% --------------------------------------------------------------

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global Nci Nco n Nri Nro et_i et_o  theta_c theta_rad theta_ini 

pl=[Nci*((ul-theta_c)/(theta_ini-theta_c))^n*(ul-1) + Nri*(1+et_i*(ul-1))*(ul^4-1^4);0];
ql =[-1;0];
pr =[Nco*((ur-theta_c)/(theta_ini-theta_c))^n*(ur-theta_c) + Nro*(1+et_o*(ur-theta_rad))*(ur^4-theta_rad^4);0];
qr=[1;0];

I am stuck with the initial boundary and need help on this. Moreover, the discontinuity part is creating problem. Any help would be appreciated.

Thanks!

-Meenal