2012-08-29 148 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

謝謝你讓我知道指針。 –

相關問題