Numerical solution of PDE with uniform initial condition
1 view (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
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'}
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
VBBV
on 2 Aug 2022
Ok. I just showed that if you change the model inputs it will give more precise results.
More Answers (0)
See Also
Categories
Find more on General PDEs 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!