第三章Euler方法的讲解.docVIP

Euler方法的讲解： 针对的方程：（一阶常微分方程处置问题） 举例： 首先： 右边函数可以写成M文件如下，但是也可以直接在命令窗口写入。 function Y=f(x,y) Y=y+cos(x); 其初值问题的解析真解命令为 Yreal=dsolve(Dy=y+cos(x),y(0)=2,x) 输出结果为： Yreal = -1/2*cos(x)+1/2*sin(x)+5/2*exp(x) 可以画出图形： x=0:0.1:4; Yreal=-0.5*cos(x)+o.5*sin(x)+2.5*exp(x); plot(x,Yreal,’-r’); legend(’Yreal’); xlabel(independent variable x); ylabel(real solution Y); 然后： 分别三种数值解法： Euler法 function y = DEEuler(f, h,a,b,y0,varvec)%这是Euler法的函数命令 %一阶常微分方程的一般表达式的右端函数：f 这里可用 f=inline，也可直接输入 % 自变量下限 a；上限 b % 函数初值 y0 % 积分步长 h % 常微分方程的变量组 varvec format long;%数据显示方式不影响计算和存储方式是指小数点后15位数字表示N = (b-a)/h; y = zeros(N+1,1); x = a:h:b; y(1) = y0; for i=2:N+1 y(i) = y(i-1)+h*Funval(f,varvec,[x(i-1), y(i-1)]); end format short; 补充讲解 Funval(f,varvec,[x(i-1), y(i-1)]) function fv = Funval(f,varvec,varval)%varvec为DEEuler中输入的变量组 var = findsym(f); %按字典顺序返回符号函数f中的所有变量名(符号) if length(var) 4 %if语句求值 if var(1) == varvec(1) fv = subs(f,varvec(1),varval(1)); else fv = subs(f,varvec(2),varval(2)); end else fv = subs(f,varvec,varval); end 其中subs(f,varvec(1),varval(1))表示符号变量代入函数的命令函数，这是matlab自带函数，比如 subs(a+b,a,4) %returns 4+b. subs(cos(a)+sin(b),{a,b},{sym(alpha),2}) %returns cos(alpha)+sin(2) 本题输入命令为 syms u v x y Y = DEEuler(-v+u+1, 0.1,0,0.5,1,[u v]) 隐式Euler法/梯形法 function y = DEimpEuler(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 % fx = Funval(f,var(1),x(i)); % gx = y(i-1)+h*fx - varvec(2); % y(i) = NewtonRoot(gx,-10,10,eps); %这里不是重点，同学们自己看 y(i)=(y(i-1)+h*cos(x(i)))/(1-h);% 隐式Euler法 y(i)=((1+h/2)*y(i-1)+0.5*h*cos(x(i-1))+0.5*h*cos(x(i)))/(1-h/2) %梯形法 end format short; 其中：NewtonRoot法编程如下 function root=NewtonRoot(f,a,b,eps) if(nargin==3) eps=1.0e-4; end f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end if(f2==0) root=b; end tol=1; fun=diff(sym(f)); fa=subs(sym(f),findsym(sym(f)),a); fb=subs(sym(f),findsym(sym(f)),b); df