(解初值问题各种方法比较2.doc
文本预览下载声明
解初值问题各种方法比较
一、实验目的：掌握了解各种解初值问题的方法，体会步长对问题解的影响。
二、实验内容：给定初值问题
其精确解为
三、实验要求：
分别按
（1）欧拉法，步长；
（2）改进的欧拉法，步长；
（3）四阶标准龙格-库塔法，步长；
编写程序求在节点处的数值解及误差，并比较各方法的优缺点。用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与上述结果进行比较。
四、实验过程：
1、（1）编写主函数。打开Editor编辑器，输入欧拉法法主程序语句：
function [h,k,X,Y,P]=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
xlabel(自变量 X), ylabel(因变量 Y)
title(用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解)
end
P=[X,Y];
以文件名Qeuler1.m保存。
（2）编写主函数。打开Editor编辑器，输入改进的欧拉法法主程序语句：
function [X,Y,n,P]= odtixing1(funfcn,x0,b,y0,h,tol)
n=fix((b-x0)/h);
X=zeros(n+1,1);
Y=zeros(n+1,1);
k=1; X(k)=x0;
Y(k,:)=y0; Y1(k,:)=y0;
% 绘图.
clc,x0,h,y0
% 产生初值.
for i=2:n+1
X(i)=x0+h;
fx0y0=feval(funfcn,x0,y0);
Y(i,:)=y0+h*fx0y0;
fxiyi=feval(funfcn,X(i),Y(i,:));
Y1(i,:)=y0+h*(fxiyi+fx0y0)/2;
% 主循环.
Wu=abs(Y1(i,:)-Y(i,:));
while Wutol
p=Y1(i,:),
fxip=feval(funfcn,X(i),p);
Y1(i,:)=y0+h*(fx0y0+fxip)/2,P1=Y1(i,:), Y(i,:)=p1;
end
x0=x0+h; y0=Y1(i,:);
Y(i,:)=y0; plot(X,Y,ro)
grid on
xlabel(自变量 X), ylabel(因变量 Y)
title(用梯形公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解)
end
X=X(1:n+1); Y=Y(1:n+1,:);
n=1:n+1; P=[n,X,Y]
以文件名odtixing1.m保存。
（3）编写主函数。打开Editor编辑器，输入四阶标准龙格-库塔法主程序语句：
function [k,X,Y,fxy,wch,wucha,P]=RK4 (funfcn,fun,x0,b,C,y0,h)
x=x0; y=y0;p=128; n=fix((b-x0)/h); fxy=zeros(p,1);
wucha=zeros(p,1); wch=zeros(p,1); X=zeros(p,1);
Y=zeros(p,length(y)); k=1; X(k)=x; Y(k,:)=y;
% 绘图.
clc,x,h,y
%计算
%fxy=fxy(:);
for k=2:n+1
x=x+h;a2=C(5); a3=C(6);
a4=C(7);b21=C(8); b31=C(9);
b32=C(10); b41=C(11);
b42=C(12); b43=C(13);c1=C(1);
c2=C(2); c3=C(3); c4=C(4);
x1=x+a2*h; x2=x+a3*h;
x3=x+a4*h;
k1=feval(funfcn,x,y);
y1=y+b21*h*k1;
k2=feval(funfcn,x1,y1);
y2=y+b31*h*k1+b32*h*k2;
k3=feval(funfcn,x2,y2);
y3=y+b41*h*k1+b42*h*k2+b43*h*k3;
k4=feval(funfcn,x3,y3);
fxy(k)=feval(fun,x);
y=y+h*
显示全部