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

12 次查看（过去 30 天）
Declan
Declan 2022-10-16
回答： Torsten 2022-10-16
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)

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2022-10-16
Ran in:
A class mate of yours seems to have the same problem:
https://de.mathworks.com/matlabcentral/answers/1827643-function-uses-the-third-order-runge-kutta-scheme?s_tid=srchtitle
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

类别

在 Help Center 和 File Exchange 中查找有关 Graphics Object Programming 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by