24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

ghw_nit

铁杆木虫 (正式写手)

应助: 14 (小学生)
金币: 7901.9
帖子: 653
在线: 175.8小时
虫号: 774851

[交流] 【求助】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

铁杆木虫 (正式写手)

应助: 14 (小学生)
金币: 7901.9
帖子: 653
在线: 175.8小时
虫号: 774851

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')