Numerical solution of PDE with uniform initial condition

1 view (last 30 days)
I have a PDE like this
With boundary and initial conditions given as
In this case, I have a uniform initial condition and when I use the ode15s solver to solve this PDE, the final result I get is always 1. So, the plot looks like the picture below but the original results are supposed to be curvy as time increases. Not sure what I am doing wrong and I hope someone can help. I have attached the code below.
clear all
clc
%%
pts = 2^10-1; tmax = 1.76;
ti = linspace(0,tmax,pts);
xi = linspace(0,1,pts);
h_init = ones(pts,1);
[t,y] = ode15s(@(t,y) example_model(t,y,xi'),ti,h_init);
tspan = [1 100 200 300 400 500];
figure,
for i = tspan
plot(xi,y(i,:));
hold on;
end
function dhdt = example_model(t,y,xi)
t
pts = length(y); h = y;
h([1 pts]) = 1; h0 = 100;
dx = xi(2)-xi(1);
% L and Lt
tau = 0.765; U0 = 0.37;
lambda = 11.6; L_cl = 0.4;
tt = t./tau; st = sqrt(tt);
term1 = -0.5.*(tt).^2;
term2 = 0.5.*sqrt(pi).*erf(st);
term3 = -st.*exp(-tt);
L = L_cl + U0.*tau.*(term1 + lambda.*(term2+term3));
Lt = -U0.*tau.*(t/tau^2 - lambda.*((exp(-tt).*st)./tau - exp(-tt)./(2.*tau.*st) + (3991211251234741.*exp(-tt))./(4503599627370496.*tau.*pi^(1/2).*st)));
if isnan(Lt)
Lt = 0;
end
% constant
C = 7.87e-7; alpha = C*h0^3/(3*L^4);
% derivatives
hx = FD_derivative(h,dx);
hxx = FD_derivative(hx,dx);
hxxx = FD_derivative(hxx,dx);
q = (h.^3).*hxxx;
qx = FD_derivative(q,dx);
dhdx = (Lt/L).*xi.*hx;
dhdt = dhdx - alpha.*qx;
end
function dydx = FD_derivative(y,dx)
N = length(y); dydx = zeros(N,1);
for i = 1:N
if i == 1
dydx(i) = -y(i) + y(i+1);
elseif i == N
dydx(i) = y(i) - y(i-1);
else
dydx(i) = -0.5*y(i-1) + 0.5*y(i+1);
end
end
dydx = dydx./dx;
end

Accepted Answer

VBBV
VBBV on 2 Aug 2022
Edited: VBBV on 2 Aug 2022
clear all
clc
%%
pts = 2^0.1; tmax = 1.76;
ti = linspace(0,tmax,50); % give suitable sample size
xi = linspace(0,1,50);
h_init = ones(50,1);
[t,y] = ode15s(@(t,y) example_model(t,y,xi'),ti,h_init);
tspan = [1 100 200 300 400 500];
figure,
color = {'b','r','g','y','m','k'}
color = 1×6 cell array
{'b'} {'r'} {'g'} {'y'} {'m'} {'k'}
for i = 1:length(tspan)
plot(xi,y(i,:),color{i});
hold on;
end
legend('t = 1','t = 100','t =200','t = 300','t = 400','t = 500')
function dhdt = example_model(t,y,xi)
t;
pts = length(y);
h = y;
h([1 pts]) = 0; % change this boundary condition here
h0 = 100;
dx = xi(2)-xi(1);
% L and Lt
tau = 0.765; U0 = 0.37;
lambda = 11.6; L_cl = 0.4;
tt = t./tau; st = sqrt(tt);
term1 = -0.5.*(tt).^2;
term2 = 0.5.*sqrt(pi).*erf(st);
term3 = -st.*exp(-tt);
L = L_cl + U0.*tau.*(term1 + lambda.*(term2+term3));
Lt = -U0.*tau.*(t/tau^2 - lambda.*((exp(-tt).*st)./tau - exp(-tt)./(2.*tau.*st) + (3991211251234741.*exp(-tt))./(4503599627370496.*tau.*pi^(1/2).*st)));
if isnan(Lt)
Lt = 0;
end
% constant
C = 7.87e-7; alpha = C*h0^3/(3*L^4);
% derivatives
hx = FD_derivative(h,dx);
hxx = FD_derivative(hx,dx);
hxxx = FD_derivative(hxx,dx);
q = (h.^3).*hxxx;
qx = FD_derivative(q,dx);
dhdx = (Lt/L).*xi.*hx;
dhdt = dhdx - alpha.*qx;
end
function dydx = FD_derivative(y,dx)
N = length(y); dydx = zeros(N,1);
for i = 1:N
if i == 1
dydx(i) = -y(i) + y(i+1);
elseif i == N
dydx(i) = y(i) - y(i-1);
else
dydx(i) = -0.5*y(i-1) + 0.5*y(i+1);
end
end
dydx = dydx./dx;
end
you can change the boundary condition to suit your model output inside the example model function. and give a reasonable step range at beginning
  2 Comments
Darsh
Darsh on 2 Aug 2022
Hi, but changing the boundary condition like you said changes the whole mathematical model for me. But I found my mistake. I forgot to introduce the third derivative boundary condition and it works better now. Thanks for the insight.
VBBV
VBBV on 2 Aug 2022
Ok. I just showed that if you change the model inputs it will give more precise results.

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!