MATLAB: I need help with this problem (finite difference method and ode)

finite differencefinite difference methodMATLABodeode45

Hello, I solved this question using ode45. But I'm having trouble solving for y(t) using finite difference method. Here is the code I have so far. Can you explain to me how to use finite difference method on this problem? Thanks.

clear all
close all
c=0.0001;
g=9.81;
%dy/dt=v
%dv/dt=-cv^3-g
%ODE system dy_dt=F(t,y)
F = @(t,y) [y(2);-c*(y(2))^3-g];
%dy-dt=f(t,y)
f = @(t,y) -y;
%initial conditions
y0 = [40;0];
%Numerical solution using ode45 
[t,yNumerical]=ode45(F,[0 5],y0);
%Analytical Solution
figure
plot(t,yNumerical(:,1),'--o','LineWidth',2)
hold on
plot(t,yNumerical(:,2),'--o','LineWidth',2)
legend('Position','Velocity')

Best Answer

How do you know that dy/dt at t=0 has to be 0 in order to reach a height of 5 m after 5 seconds ?

Hint: "ODE45" is not the correct solver for this problem, but "bvp4c".

For a finite-difference solution see

http://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/node9.html

Best wishes

Torsten.

Related Solutions

MATLAB: 4th order Runge-Kutta code that can solve for several intial conditions

Here's a version:

% Initial conditions
t0=0;
x0=0:0.2:1;
y0=1:0.2:2;
dt=0.1;
tz=1;
t_range= t0:dt:tz;
X = zeros(numel(x0),numel(t_range));
Y = zeros(numel(y0),numel(t_range));
x_rk=zeros(1,numel(t_range));
y_rk= zeros(1,numel(t_range));
for n = 1:numel(x0)
    x_rk(1)=x0(n);
    y_rk(1)=y0(n);
    for i=1:length(t_range)-1
        [x_rk(i+1), y_rk(i+1)]=rk_method(x_rk(i),y_rk(i),dt);
    end
    X(n,:) = x_rk;
    Y(n,:) = y_rk;
    %% PLots
    subplot(numel(x0)/2,2,n)
    plot (t_range, x_rk,'b-o',t_range, y_rk,'y-o','MarkerSize',5);
    grid
    title([num2str(x0(n)),'   ',num2str(y0(n))])
    legend('x','y')
end
disp('X = '), disp(X)
disp('Y = '), disp(Y)
%% Functions
function [dxdt, dydt]= f(x,y)
dxdt =  5*x-3*y;
dydt = -6*x+2*y;
end 
function [x1, y1]=rk_method(x0,y0,dt)
[k11, k12]=f(x0,y0);     
[k21, k22]=f(x0+0.5*dt*k11,y0+0.5*dt*k12); 
[k31, k32]=f(x0+0.5*dt*k21,y0+0.5*dt*k22); 
[k41, k42]=f(x0+dt*k31, y0+dt*k32); 
x1=x0+dt*((k11/6)+((k21+k31)/3)+(k41/6));
y1=y0+dt*((k12/6)+((k22+k32)/3)+(k42/6));
end

MATLAB: Time derivative of parameters within ODE solvers

function dy = ODE(t,y)
  dt = 1e-8;
  fm = myfun(t-dt);
  f = myfun(t);
  fp = myfun(t+dt);
  df = (fp - fm) / (2 * dt);
  ddf = (fp - 2 * f + fm) / dt^2;
  dy = y + f + df + ddf;
end

Best Answer

Related Solutions

MATLAB: 4th order Runge-Kutta code that can solve for several intial conditions

MATLAB: Time derivative of parameters within ODE solvers

Related Question