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