matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

10.3.2 向前欧拉方法的三种MATLAB程序

10.3.2  向前欧拉方法的三种MATLAB程序

   

(一) 向前欧拉公式的MATLAB主程序1

用向前欧拉公式(10.8)求解常微分方程初值问题(10.5)的数值解及其截断误差公式的MATLAB主程序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;

 

例10.3.2  用向前欧拉公式(10.8)求解初值问题

分别取,并将计算结果与精确解作比较,写出在每个子区间上的局部截断误差公式,画出数值解与精确解在区间上的图形.

解  输入程序

>> 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]

运行后屏幕显示分别取时,所给的初值问题在上的自变量的值构成的数组Xi  (i=1,2),利用向前欧拉公式(10.8)求出的与Xi  (i=1,2)对应的数值解Yi  (i=1, 2),Yi的绝对误差jwYi(i=1,2)和相对误差xwYi  (i=1,2)(略),每个子区间上的局部截断误差公式Ren2 =1/200*dy2, 近似解Y与精确解的图形.

 

    (二) 向前欧拉公式的MATLAB主程序2

用向前欧拉公式(10.8)求解常微分方程初值问题(10.5)的数值解及其截断误差公式的MATLAB主程序2

function [k,X,Y,P,REn]=Qeuler2(funfcn,x0,y0,b,h)

x=x0; n=fix((b-x)/h);

X=zeros(n+1,1); y=y0;

Y=zeros(n+1,1); k=1;

X(k)=x; Y(k)=y';

for k=2:n+1

X(k)=x+(k-1)*h,

fxy=feval(funfcn,x,y),

Y(k)=y+h*fxy

y=Y(k), k=k+1,plot(X,Y,'rp'),

grid,xlabel('自变量 X'), ylabel('因变量 Y')

title('用向前欧拉(Euler)公式计算dy/dx=f(x,y),y(x0)=y0在[x0,b]上的数值解')

end

k1=1:n;k=[0,k1]; P=[k',X,Y];

syms dy2,

REn=0.5*dy2*h^2;

 

例10.3.4  用向前欧拉公式(10.8)求解初值问题,取,并将计算结果与精确解作比较,写出在每个子区间上的局部截断误差公式,画出数值解与精确解在区间上的图形.

解  输入程序

>> subplot(2,1,1)

 x0=0;y0=1;b=2;n=10;h=2/10;

[k,X,Y,P,REn]=Qeuler2(@funfcn,x0,y0,b,h)

hold on

S1=1/6*(6+12*X+30*exp(2*X)).^(1/2); plot(X,S1,'b-')

title('用向前欧拉公式计算dy/dx=y-x/(3y),y(0)=1在[0,2]上的数值解')

legend('n=10时,dy/dx=y-x/(3y),y(0)=1在[0,2]上的数值解',' dy/dx =y-x/(3y),y(0)=1在[0,2]上的精确解')

hold off

jdwucY=S1-Y; jwY=S1-Y;

xwY=jwY./Y;k1=1:n;k=[0,k1];

P1=[k',X,Y,S1,jwY,xwY]

subplot(2,1,2)

n1=100;h1=2/100;

[k,X1,Y1,P1,Ren1]=Qeuler2(@funfcn,x0,y0,b,h1)

hold on

S2=1/6*(6+12*X1+30*exp(2*X1)).^(1/2); plot(X1,S2,'b-')

legend('n=100时,dy/dx=y-x/(3y),y(0)=1在[0,2]上的数值解',' dy/dx =y-x/(3y),y(0)=1在[0,2]上的精确解')

hold off

jwY1=S2-Y1;xwY1=jwY1./Y1;k1=1:n1;k=[0,k1];

P2=[k',X1,Y1,S2,jwY1,xwY1]

运行后屏幕显示取时,此初值问题在上的自变量的向量X,X1,利用向前欧拉公式(10.8)求出的与X,X1对应的数值解Y,Y1及其绝对误差向量jwY,jwY1和相对误差向量xwY,xwY1(略),每个子区间上的局部截断误差公式Ren2 =1/200*dy2, 数值解Y与精确解的图形.

 

(三) 自适应向前欧拉公式的MATLAB主程序

用自适应向前欧拉公式求解常微分方程初值问题(10.5)的数值解的MATLAB主程序2

function [H,X,Y,k,h,P]=QEuler(funfcn,x0,b,y0,tol)

%初始化.

pow=1/3;

if nargin < 5 | isempty(tol),

tol = 1.e-6;

end;

if nargin < 6 | isempty(trace),

trace = 0;

end;

x=x0; h=0.0078125*(b-x); y=y0(:);p=128; H=zeros(p,1);

X=zeros(p,1); Y=zeros(p,length(y)); k=1; X(k)=x; Y(k,:)=y';

% 绘图.

clc,x,h,y

end

% 主循环.

while (x<b)&(x+h>x)

  if x+h>b

h=b-x;

  end

   %计算斜率.

   fxy=feval(funfcn,x,y); fxy=fxy(:);

%计算误差,设定可接受误差.

delta=norm(h*fxy,'inf'); wucha=tol*max(norm(y,'inf'),1.0);

% 当误差可接受时重写解.

if delta<=wucha

x=x+h; y=y+h*fxy; k=k+1;

if k>length(X)

X=[X;zeros(p,1)]; Y=[Y;zeros(p,length(y))];

H=[H;zeros(p,1)];

      end

H(k)=h;X(k)=x;Y(k,:)=y'; plot(X,Y,'rp'), grid

xlabel('自变量 X'), ylabel('因变量 Y')

title('用向前欧拉(Euler)公式计算dy/dx=f(x,y),y(x0)=y0在[x0,b]上的数值解')

end

%更新步长.

if delta~=0.0

h=min(h*8,0.9*h*(wucha/delta)^pow);

end

end

if (x<b)

disp('Singularity likely.'), x

end

H=H(1:k); 

X=X(1:k);

Y=Y(1:k,:);

n=1:k; P=[n',H,X,Y]

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3100302.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1