应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】feval的使用
201/1返回列表
查看: 2010  |  回复: 19
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

ghw_nit

铁杆木虫 (正式写手)

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

现在在看《数值分析导论》韩渭敏译，379页中的边值问题，

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

» 猜你喜欢

» 抢金币啦！回帖就可以得到:

查看全部散金贴
1楼2011-01-07 19:38:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
f的m函数我试着写了好几个都不能通过，哪位大侠给演示一下如何使用，万分感谢
赞一下 回复此楼
2楼2011-01-07 20:04:38
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xiegangmai

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

帮助中有实例，找来看看就知道了。
赞一下 回复此楼
3楼2011-01-07 21:47:10
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
实例已经看过了，就是把方程的右边写成一个m文件，
function f=f(x,y)
f=-2*pi^2*sin(pi*x)*sin(pi*y);
调试时不能通过
赞一下 回复此楼
4楼2011-01-08 10:12:37
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xiegangmai

版主 (职业作家)
★ ★
robert2020(金币+2):辛苦了！ 2011-01-11 11:28:07
引用回帖:
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);
调试时不能通过

你是怎样调用的呢？可能是你调用的格式不对吧，错误提示是什么呢？

feval有两种用法：
[y1, y2, ...] = feval(fhandle, x1, ..., xn)； % 用函数句柄
[y1, y2, ...] = feval(function, x1, ..., xn)； % 直接用函数名

对你的函数，用函数句柄的话：
CODE:
fhandle = @f;
feval( fhandle, 0.1, 0.1 )

用函数名：
CODE:
feval( 'f', 0.1, 0.1 )
赞一下(1人) 回复此楼
5楼2011-01-08 10:28:14
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
>> 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 match and will be used
instead.
You can improve the performance of your code by using exact
name matches and we therefore recommend that you update your
usage accordingly. Alternatively, you can disable this warning using
warning('off','MATLAB:dispatcher:InexactCaseMatch').
This warning will become an error in future releases.
??? Error using ==> mtimes
Inner matrix dimensions must agree.

Error in ==> f at 2
f=-2*pi^2*sin(pi*x)*sin(pi*y);
Error in ==> poisson at 13
Fr(2:n,2:n)=h^2*feval(f,(1:n-1)*h,(1:n-1)*h);
错误的提示
赞一下 回复此楼
6楼2011-01-08 10:51:59
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xiegangmai

版主 (职业作家)
ghw_nit(金币+20):谢谢版主 2011-01-08 11:08:17
引用回帖:
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 ...

Poisson是什么东东，你自己写的函数吗？路径不对。

另外你代入的x，y分别为(1:n-1)*h，函数中f=-2*pi^2*sin(pi*x)*sin(pi*y)这样写是不对的。
CODE:
f=-2 * pi ^2 * sin( pi * x ) .* sin( pi * y );

[ Last edited by xiegangmai on 2011-1-8 at 10:59 ]
赞一下 回复此楼
7楼2011-01-08 10:57:32
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
现在的错误是这样的，
??? Subscripted assignment dimension mismatch.

Error in ==> poisson at 13
Fr(2:n,2:n)=h^2*feval(f,(1:n-1)*h,(1:n-1)*h);
poisson是一个函数
赞一下 回复此楼
8楼2011-01-08 11:04:41
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
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);
            dif=abs(temp-u(i,j));
            if(err<=diff)
                err=dif;
            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;
%surf(X,Y,U')
%xlabel('x-axis')
%ylabel('y-axis')
%zlabel('the numerical solution')
%title('Plot of the mumerical solution')
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')
s1=springtf('h=6.4f',h)
title(s1)
hold on

Err=sin(pi*X)*sin(pi*Y')-U;
subplot(1,2,2)
surf(X,Y,Err')
赞一下 回复此楼
9楼2011-01-08 11:05:13
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
这是书中的程序，求解泊松方程的程序
赞一下 回复此楼
10楼2011-01-08 11:06:14
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
我想调试一下，按着书中的程序就出现了上述的情况
赞一下 回复此楼
11楼2011-01-08 11:06:47
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
f代表泊松方程中的函数，g是边界条件
赞一下 回复此楼
12楼2011-01-08 11:09:33
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xiegangmai

版主 (职业作家)
引用回帖:
Originally posted by ghw_nit at 2011-01-08 11:09:33:
f代表泊松方程中的函数，g是边界条件

你的边界函数呢？
赞一下 回复此楼
13楼2011-01-08 11:34:09
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
边界函数是0
回复此楼
14楼2011-01-08 11:48:22
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
求解的范围是单位的正方形，边界上的函数值都是0
赞一下 回复此楼
15楼2011-01-08 11:53:22
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xiegangmai

版主 (职业作家)
★ ★ ★ ★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
robert2020(金币+5):辛苦了！ 2011-01-11 11:28:48
引用回帖:
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函数，有人可 ...

大概修改了一下，不知道结果对不对，自己验证一下。
待求解方程函数：
CODE:
function ff = f( x, y )
ff = -2 * pi ^ 2 * sin( pi * x ) * sin( pi * y' );

边界条件：
CODE:
function ff = g( x, y )
ff = 0;

Poisson函数：
CODE:
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' )

在命令窗口输入：
CODE:
ff = @f;
gg = @g;
n = 20;
u = poisson( ff, gg, n );

得到的结果：
赞一下(2人) 回复此楼
16楼2011-01-08 12:03:07
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
已经调试通过了，谢谢版主
赞一下 回复此楼
17楼2011-01-08 13:37:10
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xiegangmai

版主 (职业作家)
引用回帖:
Originally posted by ghw_nit at 2011-01-08 13:37:10:
已经调试通过了，谢谢版主

是不是误差太大了？
赞一下 回复此楼
18楼2011-01-08 14:13:18
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)
能帮我解释一下，
ff = -2 * pi ^ 2 * sin( pi * x ) * sin( pi * y' )
y为什么一定要是转置呢？
赞一下 回复此楼
19楼2011-01-08 14:32:43
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xiegangmai

版主 (职业作家)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
引用回帖:
Originally posted by ghw_nit at 2011-01-08 14:32:43:
能帮我解释一下，
ff = -2 * pi ^ 2 * sin( pi * x ) * sin( pi * y' )
y为什么一定要是转置呢？

不好意思，搞错了，正确的应该是：
CODE:
ff = -2 * pi ^ 2 * sin( pi * x' ) * sin( pi * y );

你代入的x也y分别为一个行向量，把x转置一下，再相乘得到的才是需要的矩阵。
赞一下(1人) 回复此楼
20楼2011-01-08 15:21:26
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 ghw_nit 的主题更新
201/1返回列表
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭