我正在尝试实现一个非常简单的 4 阶 Runge-Kutta 方法,用于解决 ODE y'=f(x,y)
。
我已经在 R 和 MATLAB 中实现了该算法(见下文(,但由于某种原因,在 MATLAB 中运行需要几分钟,在 R 中运行需要几毫秒。
我的问题是为什么?
似乎唯一的区别是初始化,而玩弄它似乎没有区别。
R 脚本:
# Initialise Variables ----------------------------------------------------
L = 1 #Domain of solution function
h = 0.01#step size
x0 = 0
y0 = 0
x = x0
y = y0
# Define Forcing Function -------------------------------------------------
force = function(x,y){
-16*y + 15*exp(-x)
}
# Compute Algorithm -------------------------------------------------------
for(i in 0:(L/h)){
k1 = h*force(x,y[i])
k2 = h*force(x + h/2, y[i] + k1 /2)
k3 = h*force(x + h/2, y[i] + k2 /2)
k4 = h*force(x + h , y[i] + k3 )
temp=y[i] + (1/6)*(k1 + 2*k2 + 2*k3 + k4)
y = c(y,temp)
x = x + h
i = i+1
}
t <- seq(from=0, to=L, by=h)
plot(t,y,type="l")
MATLAB 脚本:
%% Initialise Variables
y0 = 0;
x0 = 0;
L = 1; %Length of Domain of function (here y)
N = 100; %Number of steps to break the domain into
h = L/N; %Step Size
yi = y0; %intermediate value of y
xi = x0; %intermediate value of x to be ticked in algo
y = zeros(N,1); %store y values as components
x = 0:h:L; %just for plot later
%% Define Forcing Function
syms f(A,B)
f(A,B) = 15*exp(-A) - 16*B;
%% Execute Algorithm
for n = 1:1:N;
xi= h*(n-1);
k1= h*f(xi,yi);
k2= h*f(xi + h/2 , yi + k1/2);
k3= h*f(xi + h/2 , yi + k2/2);
k4= h*f(xi + h , yi + k3 );
yi= yi + (1/6)*(k1 + 2*k2 + 2*k3 + k4);
y(n,1)=yi;
end
%plot(x,y)
您遇到此问题的原因是因为您使用符号变量来执行应该是纯数字的计算。
如果按如下方式定义f
:
f = @(A,B)15*exp(-A) - 16*B;
循环几乎立即结束。还有几点说明:
- 上述语法用于定义匿名函数的函数句柄。
- 生成的
x
和y
向量具有不同的长度,因此之后您将无法plot
它们。 - 将来,您应该
profile
代码以查找性能瓶颈。
附言:
你在R中的函数定义的MATLAB等效物将非常相似:
function out = force(x)
out = 15*exp( -x(1) ) - 16*x(2);
end
或
function out = force(x,y)
out = 15*exp(-x) - 16*y;
end
。取决于输入是否作为向量出现。