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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 【求助】feval的使用
201/1返回列表
查看: 2203  |  回复: 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返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 材料/农业专业，07/08开头均可，过线就行 +3 呵唔哦豁 2026-03-23 4/200 2026-03-23 22:30 by 汪！？！
[考研] 一志愿北京化工大学 070300 学硕 336分 求调剂 +6 vv迷 2026-03-22 6/300 2026-03-23 22:00 by fuyu_
[考研] 材料292调剂 +3 橘颂思美人 2026-03-23 3/150 2026-03-23 19:36 by pswait
[考研] 07化学280分求调剂 +3 722865 2026-03-23 3/150 2026-03-23 18:48 by macy2011
[考研] 336化工调剂 +4 王大坦1 2026-03-23 5/250 2026-03-23 18:32 by allen-yin
[考研] 一志愿070300浙大化学358分，求调剂！ +4 酥酥鱼.. 2026-03-21 4/200 2026-03-23 08:12 by Iveryant
[考研] 311求调剂 +6 冬十三 2026-03-18 6/300 2026-03-22 20:18 by edmund7
[考研] 287求调剂 +8 晨昏线与星海 2026-03-19 9/450 2026-03-22 17:01 by i_cooler
[考研] 324求调剂 +6 lucky呀呀呀鸭 2026-03-20 6/300 2026-03-22 16:01 by ColorlessPI
[考研] 286分人工智能专业请求调剂愿意跨考！ +4 lemonzzn 2026-03-17 8/400 2026-03-21 22:49 by lemonzzn
[考研] 材料与化工（0856）304求B区调剂 +3 邱gl 2026-03-20 7/350 2026-03-21 19:05 by 15709483992
[考研] 求助 +5 梦里的无言 2026-03-21 6/300 2026-03-21 17:51 by 学员8dgXkO
[考研] 求调剂 +3 .m.. 2026-03-21 4/200 2026-03-21 16:25 by barlinike
[考研] 279求调剂 +5 红衣隐官 2026-03-21 5/250 2026-03-21 14:59 by lature00
[考研] 332求调剂 +3 凤凰院丁真 2026-03-20 3/150 2026-03-21 10:27 by luoyongfeng
[考研] 083200学硕321分一志愿暨南大学求调剂 +3 innocenceF 2026-03-17 3/150 2026-03-21 02:35 by JourneyLucky
[考研] 一志愿华中科技大学，080502，354分求调剂 +5 守候夕阳CF 2026-03-18 5/250 2026-03-21 01:06 by JourneyLucky
[考研] 一志愿中海洋材料工程专硕330分求调剂 +8 小材化本科 2026-03-18 8/400 2026-03-20 23:16 by JourneyLucky
[考研] 一志愿武汉理工材料工程专硕调剂 +9 Doleres 2026-03-19 9/450 2026-03-20 22:36 by JourneyLucky
[考研] 材料学求调剂 +4 Stella_Yao 2026-03-20 4/200 2026-03-20 20:28 by ms629
信息提示
关闭
请填处理意见
关闭