[理学]微分方程数值解-matlab.ppt

MATLAB ODE初值问题的数值解 PDE问题的数值解 MATLAB求常微分方程数值解的函数 基于龙格－库塔法， MATLAB求常微分方程数值解的函数，一般调用格式为： [t,y]=ode23(fname,tspan,y0) [t,y]=ode45(fname,tspan,y0) 其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为[t0,tf],表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。 例 用古典显式格式求解抛物型方程 初始条件为 边界条件为 取步长⊿x = h = 0.2 , ⊿t = k = 0.02 。 解 r = k / h2 = 0.02 / 0.22 = 0.5, 古典显式格式为 function u=gu_dian(f,a,b,c1,c2,m,n) %输入初值和U h=a/(m-1); k=b/(n-1); r=k/h^2; U=zeros(n,m); %赋边界条件 U(2:n,1)=c1; U(2:n,m)=c2; %赋初始条件 U(1,1:m)=fg(0:h:h*(m-1)); %计算内点上u的数值解U for i=2:n for j=2:(m-1) U(i,j)=(1-2*r)*U(i-1,j)+r*(U(i-1,j-1)+U(i-1,j+1)); end end % gu_dianl1.m 步长h=0.20 , k=0.02 , r = k / h2 = 0.5 a=1; b=0.20; c1=0; c2=0; m=6; n=11; U=gu_dian(fg, a,b,c1,c2,m,n) x=0:0.2:a; y=0:0.02:b; [X,Y]=meshgrid(x,y); surf(X,Y,U) % 输入U后再画图 热方程： 波动方程： -4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 -4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 -4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 -4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 -4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 -4 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 -4 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 Matrixs 椭圆方程： 特征值方程： 热方程： 波动方程： 离散方程 … ┇ ┇ ┇ ┇ … 0 , 0 0.2 , 0 0.4 , 0 1,0 0 , 0.2 ┇ 0,0.02 1 , 0.2 ┇ 1, 0.02 function y=fg(x) y=4.*(x-x.^2); (一) 矩形求解公式 由 可得 令 -----------(12) (12)式称为矩形公式(矩形法) 实际上就是Euler求解公式 (二) 梯形求解公式 由 可得 令 ------(13) 称(13)式为梯形求解公式(梯形法) 注意:(13)式是隐形公式 则梯形公式第k步的截断误差为 显然梯形法具有二阶精度 由于梯形公式为隐形公式,一般情况下不易显化 ------(14) 以上公式称为改进的Euler求解公式(改进Euler法),即 -