Info

此问题已关闭。 请重新打开它进行编辑或回答。

How to solve steady second order ODE in 1-D with discontinuos domain?

1 次查看（过去 30 天）
MEENAL SINGHAL
MEENAL SINGHAL 2020-6-19
关闭： MATLAB Answer Bot 2021-8-20
Hello all,
I am stuck with a set of steady second order differential equations in 1-D, which previously i solved using pdepe. But later i realised that the steady state equation (without time derivative) cannot be solved with PDEPE.
To add on the complexity, my equations contain discontinuos domain.
Pdes are defined as in pdes.jpg and the corresponding pdepe 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 also looked for bvp4c but didn't reach a conclusion, as to how to proceed.
Any kind of help regarding the problem is appreciated.
Thanks and regards,
Meenal

此问题已关闭。

回答（0 个）

此问题已关闭。

标签

Community Treasure Hunt

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

Start Hunting!

Translated by