📄 rk4demo.mht
字号:
From: <Saved by Windows Internet Explorer 7>
Subject:
Date: Tue, 12 May 2009 09:53:49 -0700
MIME-Version: 1.0
Content-Type: text/html;
charset="Windows-1252"
Content-Transfer-Encoding: quoted-printable
Content-Location: http://www.mece.ualberta.ca/Courses/mec390/390code/rk4demo.m
X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2900.3350
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=3DContent-Type content=3D"text/html; =
charset=3Dwindows-1252">
<META content=3D"MSHTML 6.00.6000.16825" name=3DGENERATOR></HEAD>
<BODY><PRE>% 4th order Runge-Kutta Demo
%
% diff equation: dy/dx =3D y - SIN(x) + COS(x) subject to y(0) =
=3D 1
%
% file fcn1.m function dyxy =3D fcn1(x,y)
% dydx =3D y - sin(x) + cos(x);
%
% exact solution is, y(x) =3D SIN(x) + EXP(x)
h =3D input('Enter step size: ');
printint =3D input('Enter print interval: ');
fprintf('\n\n Step size: %g print interval: %g\n\n', h, printint);
x =3D 0; % initial conditions
xmax =3D 5;
yRK =3D 1;
yeuler =3D yRK;
fprintf(' x yexact y Euler EulerError y RK4 =
RK4Error\n');
fprintf('%6.2f %10.4f\n',x,yRK);
i =3D 0;
while x < xmax
yeuler =3D yeuler + fcn1(x, yeuler) * h; % Euler solution
=20
k1 =3D fcn1(x, yRK); % 4th order Runge-Kutta =
solution
k2 =3D fcn1(x + h / 2, yRK + h * k1 / 2);
k3 =3D fcn1(x + h / 2, yRK + h * k2 / 2);
k4 =3D fcn1(x + h, yRK + h * k3);
yRK =3D yRK + (k1 + 2*k2 + 2*k3 + k4)*h / 6;
x =3D x + h;
yexact =3D sin(x) + exp(x);
Eulererror =3D (yeuler - yexact) / yexact * 100;
RKerror =3D (yRK - yexact) / yexact * 100;
i =3D i + 1;
if (rem(i,printint) =3D=3D 0) % print every printint step
fprintf('%8.3f %10.4f %12.4f %10.3f %12.4f %12.3f\n',...
x,yexact,yeuler,Eulererror,yRK,RKerror);
i =3D 0;
end;
end;
% Or try the one step Matlab solution...
[xODE45,yODE45] =3D ode45('fcn1',0,5,1);
fprintf('\n x ode45 solution error\n');
for i =3D 1:length(xODE45)
x =3D xODE45(i);
yexact =3D sin(x) + exp(x);
ODE45error =3D (yODE45(i) - yexact) / yexact * 100;
fprintf('%12.4f %12.4f %12.5f\n',x,yODE45(i),ODE45error);
end;
</PRE></BODY></HTML>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -