Clear Filters
Clear Filters

How do you make a script which uses the third- order Runge- Kutta scheme to solve for y'

18 views (last 30 days)
How do I make a function which uses
to solve
I think that the function should have the interface
function y = odesRK3(fun, t, y0)
% where y is a matrix conatining the solution,
% fun is the user-supplied function f(t,y),
% t is a row vector containing values of the independant variable t in
% ascending order,
% and y0 is a column vector of initial conditions
end
To test if the function works, im using the following code to call the function
f = @(t,y) [-2*y(1)+y(2); 2*y(1)-y(2)];
h = 0.2;
t = 0:h:1.4;
y0 = [9;3];
y = odesRK3(f, t, y0);
yExact = [4+5*exp(-3*t); 8-5*exp(-3*t)];
overallError = max(abs(y-yExact), [], 2)
plot(t, y, 'o--', t, yExact)
xlabel('t')
legend('RK3 y_1','RK3 y_2','Exact y_1','Exact y_2')
btw, im assuming that fun(t, y) returns a column vector given a scalar t and a column vector y
also, the solution at each point in the vector t must be stored in the columns of y, as in, y(: ,k) should be the solution at t(k)

Answers (1)

Torsten
Torsten on 16 Oct 2022
A class mate of yours seems to have the same problem:
Maybe you can exchange ideas.
And - together we are strong - here is the combined code:
f = @(t,y) [-2*y(1)+y(2); 2*y(1)-y(2)];
h = 0.2;
t = 0:h:1.4;
y0 = [9;3];
y = odesRK3(f, t, y0);
yExact = [4+5*exp(-3*t); 8-5*exp(-3*t)];
overallError = max(abs(y-yExact), [], 2)
overallError = 2×1
0.1760 0.1760
plot(t, y, 'o--', t, yExact)
xlabel('t')
legend('RK3 y_1','RK3 y_2','Exact y_1','Exact y_2')
function y = odesRK3(f, t, y0)
% where y is a matrix conatining the solution,
% fun is the user-supplied function f(t,y),
% t is a row vector containing values of the independant variable t in
% ascending order,
% and y0 is a column vector of initial conditions
n=length(t);
y=nan(length(y0),n);
y(:,1)=y0(:);
for k=1:n-1
h=t(k+1)-t(k);
F1=h*f(t(k),y(:,k));
F2=h*f(t(k)+h/2,y(:,k)+(F1/2));
F3=h*f(t(k)+0.75*h,y(:,k)+(0.75*F1));
y(:,k+1)=y(:,k)+(2*F1+3*F2+4*F3)/9;
end
end

Tags

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!