十常微分方程(组)求解.doc

10.1 常微分方程（组）的MATLAB符号求解 10.1.1 用MATLAB求常微分方程（组）的通解 调用格式一： S=dsolve (eqn，var) 调用格式二： S=dsolve (eqn1,eqn2, ... ,eqnm，var) 10.1.2 用MATLAB求常微分方程（组）的特解 调用格式三： S=dsolve (eqn，condition1，…，conditionn，var) 调用格式四：S=dsolve(eqn1,eqn2, ... ,eqnm，condition1，condition2，…，var) 10.1.3 线性常微分方程组的解法 求解齐次线性常微分方程组的MATLAB主程序 名为qcxxcwz.m的求解齐次线性常微分方程组的MATLAB主程序如下： function [X,E,V]=qcxxcwz(A) syms x E=eig(A); [V,n]=eig(A); X=exp(E*x)*V; 10.3 欧拉（Euler）方法的MATLAB程序 10.3.1 向前欧拉公式及其误差估计 例10.2.1 用欧拉方法求初值问题 的数值解，分别取，并计算误差，画出精确解和数值解的图形. 解 编写并保存名为Eulerli1.m的MATLAB计算和画图的主程序如下 function P=Eulerli1(x0,y0,b,h) n=(b-x0)/h; X=zeros(n,1); Y=zeros(n,1); k=1; X(k)=x0; Y(k)=y0; for k=1:n X(k+1)=X(k)+h; Y(k+1)=Y(k)+h*(X(k)-Y(k)); k=k+1; end y=X-1+2*exp(-X); plot(X,Y,mp,X,y,b-) grid xlabel(自变量 X), ylabel(因变量 Y) title(用向前欧拉公式求dy/dx=x-y，y(0)=1在[0,1]上的数值解和精确解y=x-1+2 exp(-x)) legend(h=0.075时，dy/dx=x-y，y(0)=1在[0,1]上的数值解,精确解y=x-1+2 exp(-x)) jwY=y-Y;xwY=jwY./y; k1=1:n;k=[0,k1]; P=[k,X,Y,y,jwY,xwY]; x0=0;y0=1;b=1;h=0.0750; P=Eulerli1(x0,y0,b,h) 在MATLAB工作窗口输入下面的程序 h1=0.0075; P1=Eulerli1(x0,y0,b,h1) legend(h1=0.0075时，dy/dx=x-y，y(0)=1在[0,1]上的数值解,精确解y=x-1+2 exp(-x)) 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./