Solving ODE using Euler Method and 4th Order Runge Kutta Method
87 views (last 30 days)
Show older comments
I am trying to learn how to solve differential equations provided the intial conditions, I have already made the matlab code for both the euler and runge kutta method as follows:
%Euler method
function y=elrl(t,x,n,h)
%t:Time
%x:Variables
%n:Number of variables
%h:Step
d=f(t,x);
for i = 1:n
p(i)=x(i)+h*d(i);
end
d=f(t+h,p);
for i = 1:n
q(i)=x(i)+h*d(i);
y(i)=0.5*(p(i)+q(i));
end
%Runge Kutta method
function y=rktl(t,x,n,h)
%t:Time
%x:Variables
%n:Number of variables
%h:Step
k0=f(t,x);
k1=f(t+0.5*h,x+k0*0.5*h);
k2=f(t+0.5*h,x+k1*0.5*h);
k3=f(t+h,x+k2*h);
y=x+h/6*(k0+2*k1+2*k2+k3);
My problem is that I don't know how to define the function using the initial conditions, I have the following system, can anybody help ?
9 Comments
Accepted Answer
James Tursa
on 21 Dec 2021
Edited: James Tursa
on 21 Dec 2021
You pretty much have the Euler scheme worked out, so I will help you with a vector formulation. Take this code:
for i=1:n
y0=y0+dt*y1;
y1=y1+dt*(-y0); % but as noted this line isn't quite correct
y2=y2+dt*(-y2);
t=t+dt;
p(i,:)=[y0 y1 y2];
end
Here you are filling p(i.:) at each step with the calculated y0, y1, and y2 variables which you calculate individually. However, you could code this directly as a vector like this (note that I have switched the indexes so that the state is a column vector):
for i=1:n-1
p(:,i+1) = p(:,i) + dt * f(t,p(:,i));
t=t+dt;
end
In this case the p(:,i) is a 3-element state vector of [y0;y1;y2] at the current time t, and p(:,i+1) is the 3-element state vector at the next time t+dt. It follows your scheme where p(1,i) is y0, p(2,i) is y1, and p(3,i) is y2. The only thing missing is the vector derivative function f, which could be simply this:
f = @(t,y) [y(2);-y(1);-y(3)]; % where y(1)=y0, y(2)=y1, y(3)=y2
Why did I switch the indexing? By having the states in columns, your derivative function will match what the MATLAB supplied ode functions such as ode45 expect, and it will be easy for you to double check your results by calling ode45 using the same f function. Also, it will be easier to take this vector formulation and extend it to the Modified Euler method and the RK4 scheme.
More Answers (1)
Torsten
on 21 Dec 2021
Edited: Torsten
on 23 Dec 2021
function main
tstart = 0.0;
tend = 10.0;
h = 0.01;
T = (tstart:h:tend).';
Y0 = [-1 0 1];
f = @(t,y) [y(2) -y(1) -y(3)];
Y = euler_explicit(f,T,Y0); % Euler explicit solution
plot(T,Y)
hold on
Y = runge_kutta_RK4(f,T,Y0);
plot(T,Y) % Runge Kutta solution
%hold on
%plot(T,[-cos(T) sin(T) exp(-T)]) % Analytical solution
end
function Y = euler_explicit(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for i = 2:N
t = T(i-1);
y = Y(i-1,:);
h = T(i) - T(i-1);
Y(i,:) = y + h*f(t,y);
end
end
function Y = runge_kutta_RK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for i = 2:N
t = T(i-1);
y = Y(i-1,:);
h = T(i) - T(i-1);
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
You must improve your programming skills.
So I hope you do not only copy the code and submit it as your homework, but also try to understand what's going on here.
2 Comments
Torsten
on 21 Dec 2021
The main principle you have to follow is that the variables you use in a certain line of a function must be defined or calculated in the lines before.
So you should first consider what you want the function to supply as output, what variables are known (input to the function) and what variables you need to calculate from these known quantities to produce the desired output.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!