向前欧拉公式在Matlab解微分方程初值解的问题0
fuqilin1202013.07.04浏览527次分享举报
用向前欧拉公式(10.8)求解初值问题，dy/dx=-3x+8x-7,y(0)=1,分别取n=10,n=100，并将计算结果与精确解作比较，写出在每个子区间[xk,xk+1]上的局部截断误差公式，画出数值解与精确解在区间[0,1]上的图形.
主程序：
function [h,k,X,Y,P,REn]=Qeuler1(funfcn,x0,y0,b,n,tol)
x=x0; h=(b-x)/n; X=zeros(n,1); y=y0;
Y=zeros(n,1); k=1; X(k)=x; Y(k)=y';
for k=2:n+1
fxy=feval(funfcn,x,y);
delta=norm(h*fxy,'inf');
wucha=tol*max(norm(y,'inf'),1.0);
if delta>=wucha
x=x+h; y=y+h*fxy; X(k)=x;Y(k)=y';
end
plot(X,Y,'rp')
grid
label('自变量 X'), ylabel('因变量 Y')
title('用向前欧拉(Euler)公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')
end
P=[X,Y];
syms dy2,
REn=0.5*dy2*h^2;
COMMAND WINDOW 输入:
subplot(2,1,1)
x0=0;y0=1;b=1-1.e-4;
n=100;tol=1.e-4;
[h1,k1,x1,Y1,P1,Ren1]=QEuler1(@funfcn,x0,y0,b,n,tol)
hold on
S1= 8/3*x1-29/9+38/9*exp(-3*x1), plot(x1,S1,'b-')
title('用向前欧拉公式计算dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解')
legend('n=100时，dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解',' dy/dx=8x-3y-7，y(0)=1在[0,1]上的精确解')
hold off
jdwuc1=S1-Y1; jwY1=S1-Y1;
xwY1=jwY1./S1;k1=1:n;k=[0,k1];
P1=[k',x1,Y1,S1,jwY1,xwY1]
subplot(2,1,2)
n1=10; [h2,k2,x2,Y2,P2,Ren2]=QEuler1(@funfcn,x0,y0,b,n1,tol)
hold on
S1 = 8/3*x2-29/9+38/9*exp(-3*x2), plot(x2,S1,'b-')
legend('n=10时，dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解',' dy/dx=8x-3y-7，y(0)=1在[0,1]上的精确解')
hold off
jwY2=S1-Y2;xwY2=jwY2./S1;k1=1:n1;k=[0,k1]; P2=[k',x2,Y2,S1,jwY2,xwY2]
结果出现问题：??? Error using ==> feval
Undefined function 'funfcn'.
Error in ==> C:\MATLAB6p5\work\Qeuler1.m
On line 5 ==> fxy=feval(z,x,y);
怎么改？？？
