2012-08-29 237 views
2

我有一个简单的一阶微分方程系统需要在Matlab中求解。它看起来像如下:Fortran |将函数作为其他函数的参数传递

function passing 
y0 = [0; 0]; 
t = 0:0.05:1; 
y = myprocedure(@myfunc,0.05,t,y0); 
myans = y' 
end 

function f = myfunc(t,y) 
f = [-y(1) + y(2); 
    -0.8*t + y(2)]; 
end 

function y=myprocedure(f,h,t,y0) 
n = length(t); 
y = zeros(length(y0),n); 
y(:,1) = y0; 
for i=1:n-1 
    k1 = f(t(i),y(:,i)); 
    k2 = f(t(i)+h/2, y(:,i)+h/2*k1); 
    k3 = f(t(i)+h/2, y(:,i)+h/2*k2); 
    k4 = f(t(i)+h, y(:,i)+h*k3); 
    y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4); 
end 
end 

因此,运行passing.m文件后,我有结果:

>> passing 

myans = 

    0   0 
-0.0000 -0.0010 
-0.0001 -0.0041 
-0.0005 -0.0095 
-0.0011 -0.0171 
-0.0021 -0.0272 
-0.0036 -0.0399 
-0.0058 -0.0553 
-0.0086 -0.0735 
-0.-0.0946 
-0.0169 -0.1190 
-0.0225 -0.1466 
-0.0293 -0.1777 
-0.0374 -0.2124 
-0.0469 -0.2510 
-0.0579 -0.2936 
-0.0705 -0.3404 
-0.0849 -0.3917 
-0.1012 -0.4477 
-0.1196 -0.5086 
-0.1402 -0.5746 

现在,我想的passing.m代码转换为Fortran 90的。我尝试了很多东西,但我仍然迷失了方向 - 特别是在进行数组初始化并将函数传递给其他函数时。

任何帮助表示赞赏。谢谢!

==编辑:

我成功复制上面的Matlab代码使用Fortran 90.我仍然不知道如何构建指针,但“接口”语句似乎工作。如果您查看代码并提供一些输入,我将不胜感激。谢谢。

module odemodule 
implicit none 
contains 
function odef(a,b) 
implicit none 
real::a 
real::b(2) 
real::odef(2) 
odef(1) = -b(1) + b(2) 
odef(2) = -0.8*a + b(2) 
end function odef 
subroutine rk4(f,h,t,y0,y) 
implicit none 
real,dimension(2,21)::y 
real,dimension(2)::k1,k2,k3,k4 
real::h 
real::t(21) 
real::y0(2) 
integer::i 
interface 
    function f(a,b) 
    real::a 
    real::b(2), f(2) 
    end function f 
end interface 
y(1,1)=y0(1) 
y(2,1)=y0(2) 
do i=1,20  
k1 = f(t(i),y(:,i)) 
k2 = f(t(i)+h/2, y(:,i)+h/2*k1) 
k3 = f(t(i)+h/2, y(:,i)+h/2*k2) 
k4 = f(t(i)+h, y(:,i)+h*k3) 
y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4) 
end do 
end subroutine rk4 
end module odemodule 

program passing 
    use odemodule 
implicit none 
real,parameter::dt=0.05 
real::yinit(2) 
real::y(2,21) 
integer::i 
real::t(21) 
t(1)=0.0 
do i=2,21 
t(i)=t(i-1)+dt 
end do 
yinit(1)=0 
yinit(2)=0 
call rk4(odef,dt,t,yinit,y) 
do i=1,21 
    write(*,100) y(1,i), y(2,i) 
enddo 
100 format(f7.4,3x,f7.4) 
end program passing 
+0

@HighPerformanceMark - 没问题。我花了很长时间才弄清楚了。 – mgilson

回答

3

Fortran中选择传递给其他过程(子例程或函数)的函数的最佳方法是使用过程指针。这样你就可以用通用的方式编写过程并用特定的函数调用它。这里有一些例子:How to alias a function name in FortranFunction pointer arrays in Fortran

有很多方法来初始化arrrays。例如,您可以初始化所有元素为相同的值,例如,

array = 0.0 

http://en.wikipedia.org/wiki/Fortran_95_language_features

+0

谢谢你让我知道指针。 –

相关问题