Writing a solver for the implicit Runge-Kutta method (order 4)

37 次查看（过去 30 天）
Stuart-James Burney
评论： lenin becerra merchan 2020-10-21
Hi.
I am in the process of writing a solver for the implicit Runge-Kutta (IRK) method of order 4, which will be used to solve a system of equations.
I have chosen not to include the summation inside the loop, as v=2 and it was simpler to write out the equations. Each equation for Xi is defined implicit, that is x_1 and x_2 are implicit equations. I'm unsure of how to complete my function such that it solves the implicit equations. I have thought of using fsolve, but i'm unsure of how to correctly place this command inside my loop.
Any suggestions are kindly appreciated.
Here, the values of y0 and f are provided in the file errors.m
function [tout,yout] = IRK4Solver(f,t,y0)
t = t(:); % ensure that t is a column vector
N = length(t);
h = (t(end)-t(1))/(N-1); % calculate h by assuming t gridpoints are equidistant
d = length(y0); % defines the size of y0
y0 = y0(:); % ensures that y0 is a column vector
y = zeros(d,N);% allocate the output array y
y(:,1) = y0; % assign y0 to the first column of y
c = [1/2-sqrt(3)/6;
1/2+sqrt(3)/6];
A = [1/4, 1/4-sqrt(3)/6;
1/4+sqrt(3)/6, 1/4];
b = [1/2;1/2];
%calculate the loop
for n=1:N
eq=@(xi) [xi(1:3)-(y(:,n)+h.*A(1,1).*f(t(n)+c(1).*h,xi(1:3))+h.*A(1,2).*f(t(n)+c(2).*h,xi(4:6)));
xi(4:6)-(y(:,n)+h.*A(2,1).*f(t(n)+c(1).*h,xi(1:3))+h.*A(2,2).*f(t(n)+c(2).*h,xi(4:6)))];
xistar=fsolve(eq,[1 1 1;1 1 1]);
y(:,n+1) = y(:,n)+h.*b(1).*f(t(n)+c(1).*h,xistar(1:3))+h.*b(2).*f(t(n)+c(2).*h,xistar(4:6));
end
tout = t;
yout = y;
  6 个评论
显示 4更早的评论
Walter Roberson
Walter Roberson 2020-4-26
xi_1 = y(:,n)+h.*A(1,1).*f(t(n)+c(1).*h,xi_1)+h.*A(1,2)f(t(n)+c(2).*h,xi_2);
You are missing an operation between A(1,2) and the f that follows it. You have similar problems on the next lines.
Stuart-James Burney
Stuart-James Burney 2020-4-26
编辑：Stuart-James Burney 2020-4-26
Thank you a lot, I believe I have added the missing operations. Are you able to help me complete my function by solving xi_1 and xi_2 implicitly at the beginning of my loop, so that the values can be used in the calculation of y(:,n+1)?
I have made further attempts to do this, but I have not been successful. I can include my attempts in my original post if you would like.

请先登录，再进行评论。

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

回答（1 个）

Walter Roberson
Walter Roberson 2020-4-26
Your test.m has an extra 'e' character as the very first thing.
Your lorenz is missing a multiplication.
After fixing that and adding the multiplications I mentioned above, then you reach that problem that xi_1 and xi_2 are defined implicitly.
f(t(n)+c(1).*h,xi_1)
Your f is lorenz() and the xi_1 is being passed as the x parameter, but your lorenz function assumes that its x has three elements. Therefore your xi_1 would have to be a vector of threee elements, and likewise your xi_2 would have to be a vector of three elements, for a total of six elements. But you only have two equations in those six unknowns, which is not enough to pin those down.
If you rethink your xi_1 and xi_2 as being three elements each, then you can combine the two equations, if you rethink them in terms of x = f(something,x) implies f(something,x) - x = 0 . But you will need to invent an initial guess.
After that you are left with the problem that your implicit function does not define tout or yout.
And once that is all dealt with... your other solvers invoked by errors are not defined.
  11 个评论
显示 9更早的评论
Walter Roberson
Walter Roberson 2020-4-26
You have to run and let it stop there, and then you can examine the variables.
Y was defined three lines further up at
[t_out, Y] = func(lorenz,t,y0);
lenin becerra merchan
lenin becerra merchan 2020-10-21
how can i do to include this expression
% Define ODE function
f = @(t, y) [ y(2) , -10*(1-y(1)^2)*y(2)-y(1) ];
y0 = 2

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by