Writing a solver for the implicit Runge-Kutta method (order 4)
18 次查看(过去 30 天)
显示 更早的评论
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 个评论
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.
回答(1 个)
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 个评论
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
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
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!