Numerical solution of PDE with uniform initial condition
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
VBBV
2022-8-2
编辑:VBBV
2022-8-2
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 个评论
VBBV
2022-8-2
Ok. I just showed that if you change the model inputs it will give more precise results.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 General PDEs 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!