Numerical Integration of a wave equation, using central differences to estimate dp/dx and using RK4 to integrate WRT time.

10 views (last 30 days)
PLEASE NOTE THAT THE FOLLOWING POST IS LONG, IF YOU HAVE NO INTEREST IN MY ANALYSIS AND JUST THE CODE, FEEL FREE TO JUST READ THE CODE AND PROVIDE RESPONSES.
Hi there, I am currently doing my dissertation on simulation of flow and the propagating waves in a diesel engine manifold. My task at the moment is to simulate a single propagating wave through a pipe of one meter. I could just copy and paste the solution and ask for help and solutions. But I am more interested in having a good discussion on the topic as I am thoroughly enjoying it. But obviously I am in need of help at this time.
So here is an outline of my approach before any code is mentioned. The wave I am propagating is:
(dp/dt)= -c*(dp/dx)
where C is a constant. To begin the analysis I started of an initial condition of
P(x) = Kexp((ax)^2) , where K and a are constants to be determined.
My first question is.. do we need to know the initial condition in order to simulate? I seem to have no way of doing the calculations without it. So I just chose this function. But in practice I would have the user supply their own function.
Okay, now I applied the Boundary conditions in order to find the unknowns. I assumed the pipe was closed at the end, and I obtained a final expression for the initial condition. Okay now I created my own script to calculate numerically dp/dx , using central differences approximaions. The code is placed below, and works for my proposed function.
dx=0.01
P0=1000
x=0:dx:1
dF= zeros(0,length(x));
F= @(x) P0*exp((x)*(log(3)).^0.5)
for i=2:(length(x)-1) % calculation loop
F1(i+1) = F(x(i)+dx);
F2(i-1) = F(x(i)-dx);
dF(i) = ((F1(i+1)-F2(i-1))/(2*dx)) ; % main equation
end
plot(x(:,1:length(dF)),dF,'r')
I have tried the script with multiple functions and it works great. Apart from at the ends of the pipe, as you cant use central differences there. I will modify this later with a one sided difference at the pipe ends. But for now this will do.
Okay now we have calculated dp/dx , and we have called this dF, we now have the following equation to numerically solve:
dp/dt = -C*dF
I then have created a fourth order runge-kutte method to integrate with repect to time. This code is placed below:
dt=0.01; % step size
t = 0:dt:10; % Calculates upto y(3)
Pt = zeros(1,length(t));
Pt(1) = 1000 % initial condition
Fx = @ (x)dF(i) % change the function as you desire
for n=1:(length(t)-1) % calculation loop
k1 = Fx(Pt(n));
k2 = Fx((Pt(n)+0.5*dt*k1));
k3 = Fx((Pt(n)+0.5*dt*k2));
k4 = Fx((Pt(n)+k3*dt));
Pt(n+1) = Pt(n) + (1/6)*(k1+2*k2+2*k3+k4)*dt'; % main equation
end
plot(t(:,1:length(Pt)),Pt,'r')
% code
end
Okay, first notice that In my runge-Kutte code I entered the intial condition as 1000, which is incorrect. As the initial condition is a function of x, not a constant. When I enter this initial condition I just get a straight line on my graph of t,Pt. However when I enter in the initial condition function it just gives me an error saying:
??? The following error occurred converting from
function_handle to double:
Error using ==> double
Conversion to double from function_handle is not
possible.
Okay so does anyone have any idea where I am going wrong? I have spent 22 hours on this now and Im just not getting anywhere. Im unsure if my code is incorrect, method, or analysis. I understand this forum is based on matlab coders so many people will only be able to give me code advice and not on the actual method, But I thought it would beneficial for you to fully understand the problem. Also anyone who is interested in this area of simulation, and engineering, if you want to get in touch with me for a discussion on the topic Im happy to do so.
Thank you

Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!