Numerical solution of PDE with uniform initial condition
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
VBBV
le 2 Août 2022
Modifié(e) : VBBV
le 2 Août 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 commentaires
VBBV
le 2 Août 2022
Ok. I just showed that if you change the model inputs it will give more precise results.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur General PDEs dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!