Feedback and help with my code that uses ODE15s

2 次查看(过去 30 天)
Hi all, I've developed some code that solves a moving boundary PDE with ODE15s, that I'd quite like some feedback on as it's still not producing the answers I would expect it to do from the available theory. In order to be able to use ODE15s, I have spatially discretised using the method of lines. On the fixed domain , the Original PDE system is;`
where . For the initial and boundary conditions, we have , , and . I've created two functions, one that solves for S and n using ODE15s called CellsODE, and then one that simulates called CellVelocity. Due to stability conditions, I've implemented an upwind finite difference scheme where the a forward discretisation is used for the derivative in the first equation. The other spatial second order derivatives use a central discretisation, with a forward and backward one at each boundary respectively. My velocity function is,
%CELLVELOCITY Finds V vector given n,s and parameters
function V = CellVelocity(n,s,dx)
nm=[NaN;n(1:end-1)]; %n_{i-1}
np = [n(2:end);NaN]; %n_{i+1}
dndx=[NaN;np(2:(end-1))-nm(2:(end-1));3*n(end)-4*n(end-1)+n(end-2)]/2*dx; %dndx derivative
V=[0;-dndx(2:end)]/s; %Definition of the velocity, V=0 at the boundary
end
The CellsODE function (for the first and last equation is)
%Extract n and s from the input.
s=input(space_steps+1);
n=[input(1:space_steps-1);0]; %n=0 at the boundary
x=linspace(0,1,space_steps).'; %discrete mesh
V = CellVelocity(n,s,dx); %Calling for the velocity
np = [n(2:end);NaN]; %n_{i+1} %Declaring the indexes
vm = [NaN;V(1:end-1)]; %v_{i-1}
vp = [V(2:end);NaN]; %v_{i+1}
dndx = [np(1:end-1)-n(1:end-1);n(end)-n(end-1)]/dx; %dndx derivative, first order forward
dvdx = [-V(3)+4*V(2)-3*V(1);vp(2:end-1)-vm(2:end-1);3*V(end)-4*V(end-1)+V(end-2)]/2*dx; %dvdx derivative, second order
dndt = (x*V(end)-V).*dndx/s -n.*dvdx/s +n.*(1-n); %ode for n
dsdt = V(end); %ode for s
output=[dndt;dsdt];
end
The 'Command script' that defines all the variables and calls CellsODE etc is;
T = 1000000 ; %Time
space_steps = 200; %Space steps
x = linspace(0,1,space_steps).'; %x vector
dx = 1/(space_steps-1); %space step
s0 = 1; %initial s
n0 = 0.6*(1-x.^2); %initial conditions
T_span=[0;T]; %T_vector
[T_out,Y_out] = ode15s(@(t,y) CellsODE(t,y,space_steps,dx),T_span,[n0;s0]);
My main interest is finding the function , and one can plot the adaptvie time mesh against the found by doing
Boundary = Y_out(:,end);
plot(T_out, Boundary);
Theory tells us that S for large times should be linear and the gradient should be equal to , however mine is coming out to be something very small and completely dissimilar. This is the reason why I think there is something wrong with my code, as it does not match up with the theory. One can calculate the gradient there by simply doing
grad = gradient(Boundary,T_out);
I've tried many things, I've tried changing all the discretisations to be other ones, Changing the positioning of the variables S in the code and everything under the sun. What I have noticed is that changing the number of mesh points really changes the output, which I don't think should happen. I've attached an image of the Boundary, one can see it almost goes linear but if you closer then it isn't. Thanks for reading.

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by