ch9MATLAB在数学中的应用解析.ppt

* ? for k=1:m for i=1:n s=0; for j=1:i-1 s=s+a(i,j)*y(j);end s1=0; for j=i+1:n s1=s1+a(i,j)*x(j);end y(i)=(1-w)*x(i)+w*(b(i)-s-s1)/a(i,i); end x=y; end z=x * ? 例如，给定如下方程组 给定松弛因子w=1.45，在Matlab命令窗口中调用如图9-7所示。 * ? 9.5 常微分方程的数值解 9.5.1 单步法 1.一阶常微分问题 一个常微分方程的原型可以描述为： 为了能求出常微分方程在各离散点的数值解，可以用差商近似代替初值问题的导数。 * ? 令h=x1-x0，则有 y(x1)=y(x0)+hf(x0,y0)，再利用y1 及f(x1,y1)可得到y(x2)=y(x1)+hf(x1,y1)，一般地，有 y(xn+1)=y(xn)+hf(xn,yn) ，n=0,1,…,N-1，h= xn+1-xn。 可以记y(xn+1)= yn+1，则有yn+1=yn+hf(xn,yn)。 * ? 2.Euler方法 刚才讨论的公式就是Euler方法 yn+1=yn+hf(xn,yn)，n=0,1,2,…,N-1。 算法描述如下 function[x,y]=naeuler(dyfun,xspan,y0,h) %用途：Euler格式解常微分方程 %格式：[x,y]=naeuler(dyfun,xspan,y0,h) dyfun为函数f(x,y),xspan为求解区 %间[x0,xN],y0为初值y(x0),h为步长，x返回节点,y返回数值解 * ? x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end x=x;y=y; * ? 3.改进的Euler方法 yn+1=yn+h/2[f(xn,yn)+ f(xn+1,yn+1)]，n=0,1,2,…,N-1。 算法描述如下 function[x,y]=naeuler2(dyfun,xspan,y0,h) %用途：改进Euler格式解常微分方程 %格式：[x,y]=naeuler2(dyfun,xspan,y0,h) dyfun为函数f(x,y),xspan为求解区% 间[x0,xN],y0为初值y(x0),h为步长，x返回节点,y返回数值解 * ? x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; k2=feval(dyfun, x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end x=x;y=y; * ? 4.隐式Euler公式 yn+1=yn+h f(xn+1,yn+1)，n=0,1,2,…,N-1。 算法描述如下 function[x,y]=naeulerb(dyfun,xspan,y0,h) %用途：隐式Euler格式解常微分方程 %格式：[x,y]=naeulerb(dyfun,xspan,y0,h) dyfun为函数f(x,y),xspan为求解区% 间[x0,xN],y0为初值y(x0),h为步长，x返回节点,y返回数值解 x=xspan(1):h:xspan(2); y(1)=y0; * ? for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h); end x=x;y=y; function y=iter(dyfun,x,y,h) y0=y;e=1e-4;K=1e+4; y=y+h*feval(dyfun,x,y); y1=y+2*e;k=1; while abs(y-y1)e y1=y; y=y0+ h*feval(dyfun,x,y); k=k+1; if kK,error(迭代发散);end end * ? 5.变形Euler公式 yn+1=yn+h f(xn+h/2,yn+h/2f(xn,yn))，n=0,1,2,…,N-1 算法描述如下 * ? function[x,y]=naeuler1(dyfun,xspan,y0,h) x=xspan(1):h:xs