RK23 with given coefficients

2022 12 18
1 个回答
回答已采纳
更新时间：2022 12 18
3 次查看（30 天）

0 个投票

I have an IVP; ; ;
and
I want to use RK23, but it seems my 'C' coefficient is unused as the RHS of my equation, f(u), only has one argument i.e u.
function udot = eqns(u)
udot=(-sin(2*pi*u));
end
I know it has to be used somewhere but I can't figure this out. Below is my code so far
function [u, h] = func_particle_rk3_generic(u_0, tmax, N)
a = zeros(3,3);b = zeros(1,3);c = zeros(3,1);
h = tmax/N;
u = u_0;
t = [0,tmax];
for i = 1:N
t(i+1) = i*h;
k(:,1)=eqns(u);
k(:,2)=eqns(u+h*(a(2,1)*k(:,1)));
k(:,3)=eqns(u+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)));
u = u + h*(b(1)*k(:,1) + b(2)*k(:,2) + b(3)*k(:,3));
end
It doesn't gives me an error, but it doesn't also gives me the accuracy I want to get when compared to the uExact

3 个评论
显示 1更早的评论 隐藏 1更早的评论

Torsten
Torsten 2022-12-18
编辑：Torsten 2022-12-18
a = zeros(3,3);b = zeros(1,3);c = zeros(3,1);
??? That's not RK23.
If you are only interested in u at t=tend, your code should work. But usually, one wants to see the time development of u:
u(1) = u0;
...
u(i+1) = u(i) + h*(b(1)*k(:,1) + b(2)*k(:,2) + b(3)*k(:,3));
lord_eul3r
lord_eul3r 2022-12-18
a = zeros(3,3);b = zeros(1,3);c = zeros(3,1);
these are representations of the RK23 Butcher Tableau(see attached photo on question)
in a separate script, I'd call them out as:
a = zeros(3,3);b = zeros(1,3);c = zeros(3,1);
a(2,1) = 1/3;a(3,2) = 2/3;
b(1) = 1/4;b(3) = 3/4;
c(2) = 1/3;c(3) = 2/3;
here, my concern is the 'C' values above which i didn't get to ues or don't know when to use.
lord_eul3r
lord_eul3r 2022-12-18
I think I understand why I don't need the C coefficient (IN THIS PARTICULAR CASE).
Because f doesn't depend explicitly on t, only on u. I stands corrected though

请先登录，再进行评论。

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

 采纳的回答

Torsten
Torsten 2022-12-18
Ran in:

0 个投票

Ran in:
Yes, it's unused in your equation.
But you could generalize your code as to include it:
u0 = 0.45;
f = @eqns;
f_exact = @(t)1/pi*atan(exp(-2*pi*t)*tan(pi*u0));
tmax = 0.5;
N = 100;
[T,U] = func_particle_rk3_generic(f,u0,tmax,N);
hold on
plot(T,U)
plot(T,f_exact(T))
grid on
function [t,u] = func_particle_rk3_generic(f,u0,tmax,N)
a = zeros(3,3);
b = zeros(1,3);
c = zeros(3,1);
a(2,1) = 1/3;
a(3,2) = 2/3;
b(1) = 1/4;
b(3) = 3/4;
c(2) = 1/3;
c(3) = 2/3;
h = tmax/N;
u(1) = u0;
t(1) = 0.0;
for i = 1:N-1
k(:,1) = f(t(i),u(i));
k(:,2) = f(t(i)+c(2)*h,u(i)+h*(a(2,1)*k(:,1)));
k(:,3) = f(t(i)+c(3)*h,u(i)+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)));
u(i+1) = u(i) + h*(b(1)*k(:,1) + b(2)*k(:,2) + b(3)*k(:,3));
t(i+1) = t(i) + h;
end
end
function udot = eqns(t,u)
udot=(-sin(2*pi*u));
end

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Dates and Time 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by