资源描述:
《MATLAB语言程序设计基础.docx》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库。
1、%第六章微分方程问题的解法%微分方程的解析解方法%常微分方程问题的数值解法%微分方程问题算法概述%四阶定步长Runge-Kutta算法及MATLAB实现%一阶微分方程组的数值解%微分方程转换%特殊微分方程的数值解%边值问题的计算机求解%偏微分方程的解%6.1微分方程的解析解方法%y=dsolve(f1,f2,…,fm,'x')%symst;u=exp(-5*t)*cos(2*t+1)+5;%uu=5*diff(u,t,2)+4*diff(u,t)+2*u%symsty;%y=dsolve(['D4y+1
2、0*D3y+35*D2y+50*Dy+24*y='...%'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'],'y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')%[x,y]=dsolve('D2x+2*Dx=x+2*y-exp(-t)',...%'Dy=4*x+3*y+4*exp(-t)')%symstx;%x=dsolve('Dx=x*(1-x^2)+1')%Warning:Explicitsolutioncoul
3、dnotbefound;implicitsolutionreturned.%>InD:MATLAB6p5toolboxsymbolicdsolve.matline292%x=%t-Int(1/(a-a^3+1),a=``..x)+C1=0%故只有部分非线性微分方程有解析解。%6.2微分方程问题的数值解法%6.2.1微方程问题算法概述%Euler算法%%function[outx,outy]=MyEuler(fun,x0,xt,y0,PointNum)%fun表示f(x,y);x0,xt:自变量的
4、初值和终值;%y0:函数在x0处的值,其可以为向量形式;%PointNum表示自变量在[x0,xt]上取的点数%ifnargin<5
5、PointNum<=0%PointNum默认值为100%PointNum=100;%end%ifnargin<4%y0默认值为0%y0=0;%end%h=(xt-x0)/PointNum;%计算步长h%x=x0+[0:PointNum]'*h;%自变量数组%y(1,:)=y0(:)';%将输入存为行向量,输入为列向量形式%fork=1:PointNum%f=feval(f
6、un,x(k),y(k,:));f=f(:)';%计算f(x,y)在每个迭代点的值%y(k+1,:)=y(k,:)+h*f;%对于所取的点x迭代计算y值%end%outy=y;outx=x;%plot(x,y)%画出方程解的函数图%end%[x1,y1]=MyEuler('myfun',0,2*pi,1,16);%myfun=incline'Dy=sin(x)+y'%function[Xout,Yout]=MyEulerPro(fun,x0,xt,y0,PointNumber)%%MyEulerPro用
7、改进的欧拉法解微分方程%ifnargin<5
8、PointNumber<=0%%PointNumer默认值为100%PointNumer=100;%end%ifnargin<4%y0默认值为0%y0=0;%end%h=(xt-x0)/PointNumber;%%计算所取的两离%散点之间的距离%x=x0+[0:PointNumber]'*h;%%表示出离散的自变量x%y(1,:)=y0(:)';%fori=1:PointNumber%迭代计算过程%f1=h*feval(fun,x(i),y(i,:));f1
9、=f1(:)';%f2=h*feval(fun,x(i+1),y(i,:)+f1);f2=f2(:)';%y(i+1,:)=y(i,:)+1/2*(f1+f2);%end%Xout=x;Yout=y;%6.2.2四阶定步长Runge-Kutta算法及MATLAB实现%6.2.3一阶微分方程组的数值解%6.2.3.1四阶五级Runge-Kutta-Felhberg算法%6.2.3.2基于MATLAB的微分方程%求解函数%格式1:直接求解%[t,x]=ode45(Fun,[t0,tf],x0)%格式2:带有
10、控制参数%[t,x]=ode45(Fun,[t0,tf],x0,options)%格式3:带有附加参数%[t,x]=ode45(Fun,[t0,tf],x0,options,p1,p2,…)%[t0,tf]求解区间,x0初值问题的初始状态变量。%描述需要求解的微分方程组:%不需附加变量的格式%functionxd=funname(t,x)%可以使用附加变量%functionxd=funname(t,x,flag,p1,p2,…)%t是时间变