Solving ODE using Euler Method and 4th Order Runge Kutta Method

87 views (last 30 days)
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
Adnane Ait Zidane
Adnane Ait Zidane on 21 Dec 2021
How is possible to define the function in a different script and use the euler method in other another and combine them? I have just started using matlab so i am not very proefficient.
Thank you.
Adnane Ait Zidane
Adnane Ait Zidane on 21 Dec 2021
Am I understanding correctly?
I have to save the initial value of y0 to use it to calculate y1 ?

Sign in to comment.

Accepted Answer

James Tursa
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.
  1 Comment
Adnane Ait Zidane
Adnane Ait Zidane on 21 Dec 2021
I have done the following, I get an error though. Both are 1x1 doubles so where is the problem coming from?
%Euler Method
function y=elrl(y,h)
n=1000; %Number of iterations
t=0; %Initial value
dt=0.01; %Step size
y0=-1;y1=0;y2=1;
for i=1:n
y0save=y0;
y0=y0+dt*y1;
y1=y1+dt*(-y0save);
y2=y2+dt*(-y2)
p(i,:)=[y0 y1 y2];
end
for i=1:n-1
f = @(t,y) [y(2);-y(1);-y(3)];
p(:,i+1)=p(:,i)+dt*f(t,p(:,i));
t=t+dt;
end
figure
plot (transpose(0.01:0.01:10),p)

Sign in to comment.

More Answers (1)

Torsten
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
Adnane Ait Zidane
Adnane Ait Zidane on 21 Dec 2021
I would like to understand the logic behind putting the expression of the function into that order, it would be helpful if you can explain.
Thank you.
Torsten
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.

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!