24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 2245  |  回复: 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 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 298求调剂 +4 种圣赐 2026-03-28 4/200 2026-03-29 08:42 by q1092522407
[考研] 305求调剂 +8 RuiFairyrui 2026-03-28 8/400 2026-03-29 08:22 by fmesaito
[考研] 一志愿211 初试270分 求调剂 +8 谷雨上岸 2026-03-23 9/450 2026-03-29 01:00 by 我是小康
[考研] 一志愿太原理工安全工程300分,求调剂 +5 0857求调剂. 2026-03-24 6/300 2026-03-28 22:04 by zhq0425
[考研] 化学0703 调剂 306分 一志愿211 +4 26要上岸 2026-03-28 4/200 2026-03-28 15:30 by 1018329917
[考研] 一志愿北京工业大学,324分求调剂 +4 零八# 2026-03-28 4/200 2026-03-28 15:01 by 17865157980
[考研] 304求调剂 +6 曼殊2266 2026-03-27 6/300 2026-03-28 14:10 by 唐沐儿
[考研] 0703化学求调剂 +9 奶油草莓. 2026-03-22 10/500 2026-03-28 13:30 by 唐沐儿
[考研] 266分,求材料冶金能源化工等调剂 +7 哇呼哼呼哼 2026-03-27 9/450 2026-03-28 12:22 by zllcz
[考研] 材料与化工考研调剂 +17 孅華 2026-03-22 17/850 2026-03-28 08:35 by WYUMater
[考研] 07化学280分求调剂 +10 722865 2026-03-23 10/500 2026-03-27 15:51 by Plutoqq
[考研] 08开头275求调剂 +4 拉谁不重要 2026-03-26 4/200 2026-03-27 14:12 by Delta2012
[考研] 316求调剂 +5 Pigcasso 2026-03-24 5/250 2026-03-27 12:10 by zhshch
[考研] 总分322求生物学/生化与分子/生物信息学相关调剂 +5 星沉uu 2026-03-26 6/300 2026-03-26 19:02 by macy2011
[考研] 中国科学院深圳先进技术研究院-光纤传感课题组招生-中国科学院大学、深圳理工大学联培 +5 YangTyu1 2026-03-26 5/250 2026-03-26 18:27 by 猫咪猫咪呀
[考研] 一志愿哈工大,085400,320,求调剂 +4 gdlf9999 2026-03-24 4/200 2026-03-25 23:01 by boxking200
[考研] 085404电子信息284分求调剂 +4 13659058978 2026-03-24 4/200 2026-03-24 12:15 by syl20081243
[考博] 26申博自荐 +3 whh869393 2026-03-24 3/150 2026-03-24 09:55 by 21018060
[考研] 一志愿国科过程所081700,274求调剂 +3 三水研0水立方 2026-03-23 3/150 2026-03-23 23:11 by MajorWen
[考研] 280分求调剂 一志愿085802 +4 PUMPT 2026-03-22 7/350 2026-03-22 22:13 by 星空星月
信息提示
请填处理意见