[交流] 【求助】feval的使用

Originally posted by ghw_nit at 2011-01-07 20:04:38:
f的m函数我试着写了好几个都不能通过，哪位大侠给演示一下如何使用，万分感谢

Originally posted by ghw_nit at 2011-01-08 10:12:37:
实例已经看过了，就是把方程的右边写成一个m文件，
function f=f(x,y)
f=-2*pi^2*sin(pi*x)*sin(pi*y);
调试时不能通过

fhandle = @f;
feval( fhandle, 0.1, 0.1 )

feval( 'f', 0.1, 0.1 )

Originally posted by ghw_nit at 2011-01-08 10:51:59:
>> Poisson('f','g',16,1e-8,10000)
Warning: Could not find an exact (case-sensitive) match for 'Poisson'.
C:\Documents and Settings\lenovo\My Documents\MATLAB\poisson.m is a case-insensitive ...

f=-2 * pi ^2 * sin( pi * x ) .* sin( pi * y );

Originally posted by ghw_nit at 2011-01-08 11:09:33:
f代表泊松方程中的函数，g是边界条件

Originally posted by ghw_nit at 2011-01-07 19:38:24:
现在在看《数值分析导论》韩渭敏译，379页中的边值问题，

在对这个偏微分方程的求解过程中用到一条语句
feval(f,(1:n-1)*h,(1:n-1)*h)
注释中说明f方程的右边的函数值，我不知道如何写这个f的m函数，有人可 ...

function ff = f( x, y )
ff = -2 * pi ^ 2 * sin( pi * x ) * sin( pi * y' );

function ff = g( x, y )
ff = 0;

function U = poisson( f, g, n, tol, max_it )
if nargin < 5
    max_it = 10000;
end
if nargin < 4
    tol = 1e-5;
end

n1 = n + 1;
h = 1 / n;
toln = ( h ^ 2 ) * tol;
Fr = zeros( n, n );
Fr( 2 : n, 2 : n ) = h ^ 2 * feval( f, ( 1 : n - 1 ) * h, ( 1 : n - 1 ) * h );

U = zeros( n1, n1 );
U( 1, 1 : n1 ) = feval( g, 0, ( 0 : n ) * h );
U( n1, 1 : n1 ) = feval( g, 1, ( 0 : n ) * h );
U( 1 : n1, 1 ) = feval( g, ( 0 : n ) * h, 0 );
U( 1 : n1, n1 ) = feval( g, ( 0 : n ) * h, 1 );

rel_err = 1;
itnum = 0;
while( ( rel_err > toln ) && ( itnum <= max_it ) )
    err = 0;
    umax = 0;
    for j = 2 : n
        for i = 2 : n
            temp = ( U( i + 1, j ) + U( i - 1, j ) + U( i, j + 1 ) + U( i, j - 1 ) ) / 4 - Fr( i, j );
            diff = abs( temp - U( i, j ) );
            if ( err <= diff )
                err = diff;
            end
            U( i, j ) = temp;
            temp = abs( temp );
            if( umax <= temp )
                umax = temp;
            end
        end
    end
    itnum = itnum + 1;
    rel_err = err / umax;
end

X = ( 0 : h : n * h )';
Y = X;
subplot( 1, 2, 1 )
surf( X, Y, U' )
xlabel( 'x-axis' )
ylabel( 'y-axis' )
zlabel( 'The numerical solution' )

Err = sin( pi * X ) * sin( pi * Y' ) - U;
subplot( 1, 2, 2 )
surf( X, Y, Err' )

ff = @f;
gg = @g;
n = 20;
u = poisson( ff, gg, n );

Originally posted by ghw_nit at 2011-01-08 13:37:10:
已经调试通过了，谢谢版主

Originally posted by ghw_nit at 2011-01-08 14:32:43:
能帮我解释一下，
ff = -2 * pi ^ 2 * sin( pi * x ) * sin( pi * y' )
y为什么一定要是转置呢？

ff = -2 * pi ^ 2 * sin( pi * x' ) * sin( pi * y );