向前歐拉公式在Matlab解微分方程初值解的問題0

fuqilin1202013.07.04瀏覽527次分享舉報(bào)

用向前歐拉公式(10.8)求解初值問題，dy/dx=-3x+8x-7,y(0)=1,分別取n=10,n=100，并將計(jì)算結(jié)果與精確解作比較，寫出在每個(gè)子區(qū)間[xk,xk+1]上的局部截?cái)嗾`差公式，畫出數(shù)值解與精確解在區(qū)間[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)公式計(jì)算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的數(shù)值解')

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('用向前歐拉公式計(jì)算dy/dx=8x-3y-7，y(0)=1在[0,1]上的數(shù)值解')

legend('n=100時(shí)，dy/dx=8x-3y-7，y(0)=1在[0,1]上的數(shù)值解',' 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時(shí)，dy/dx=8x-3y-7，y(0)=1在[0,1]上的數(shù)值解',' 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]

結(jié)果出現(xiàn)問題：??? Error using ==> feval

Undefined function 'funfcn'.

Error in ==> C:\MATLAB6p5\work\Qeuler1.m

On line 5 ==> fxy=feval(z,x,y);

怎么改？？？

總結(jié)

以上是生活随笔為你收集整理的matlab 向前欧拉公式,向前欧拉公式在Matlab解微分方程初值解的问题的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。