24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 2251  |  回复: 19

ghw_nit

铁杆木虫 (正式写手)


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

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

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

» 猜你喜欢

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

查看全部散金贴

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)


f的m函数我试着写了好几个都不能通过,哪位大侠给演示一下如何使用,万分感谢
2楼2011-01-07 20:04:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
引用回帖:
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的回帖
★ ★
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 )

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的回帖
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的回帖
引用回帖:
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的回帖
★ ★ ★ ★ ★ ★
小木虫(金币+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 );

得到的结果:
16楼2011-01-08 12:03:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ghw_nit

铁杆木虫 (正式写手)


已经调试通过了,谢谢版主
17楼2011-01-08 13:37:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
引用回帖:
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的回帖

小木虫(金币+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转置一下,再相乘得到的才是需要的矩阵。
20楼2011-01-08 15:21:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 ghw_nit 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 材料化工340求调剂 +3 jhx777 2026-03-30 3/150 2026-03-30 17:54 by JourneyLucky
[考研] 318求调剂 +7 陈晨79 2026-03-30 7/350 2026-03-30 10:49 by 探123
[考研] 311求调剂 +10 lin0039 2026-03-26 10/500 2026-03-30 10:26 by herarysara
[考研] 一志愿:西北大学,英一数一408-284分求调剂 +4 12.27 2026-03-27 4/200 2026-03-29 14:40 by zhshch
[考研] 329求调剂 +10 钮恩雪 2026-03-25 10/500 2026-03-29 13:32 by peike
[考研] 340求调剂 +6 Amber00 2026-03-26 6/300 2026-03-29 12:06 by 无际的草原
[考研] 0856求调剂 +13 zhn03 2026-03-25 14/700 2026-03-29 08:13 by fmesaito
[考研] 330分求调剂 +5 qzenlc 2026-03-29 5/250 2026-03-29 07:37 by 无际的草原
[考研] 学硕274求调剂 +9 Li李鱼 2026-03-26 9/450 2026-03-28 21:42 by bymhappy
[考研] 312,生物学求调剂 +3 小译同学abc 2026-03-28 3/150 2026-03-28 15:32 by 落睿可思
[考研] 0856,材料与化工321分求调剂 +12 大馋小子 2026-03-27 13/650 2026-03-28 10:56 by self2008
[有机交流] 高温高压反应求助 10+4 chibby 2026-03-25 4/200 2026-03-27 21:08 by BT20230424
[考研] 285求调剂 +4 AZMK 2026-03-27 7/350 2026-03-27 20:59 by AZMK
[考研] 336材料求调剂 +7 陈滢莹 2026-03-26 9/450 2026-03-27 00:20 by wxiongid
[考研] 327求调剂 +7 prayer13 2026-03-23 7/350 2026-03-26 20:48 by 不吃魚的貓
[考研] 一志愿河工大 081700 276求调剂 +4 地球绕着太阳转 2026-03-23 4/200 2026-03-26 14:27 by zzll406
[考研] 309求调剂 +4 gajsj 2026-03-25 5/250 2026-03-26 00:27 by Dyhoer
[考研] 303求调剂 +6 蓝山月 2026-03-25 6/300 2026-03-25 22:47 by 418490947
[考研] 07化学303求调剂 +5 睿08 2026-03-25 5/250 2026-03-25 22:46 by 418490947
[考研] 086003食品工程求调剂 +6 淼淼111 2026-03-24 6/300 2026-03-25 10:29 by 3Strings
信息提示
请填处理意见