【MATLAB】实验五：数值微积分与方程数值求解

实验五数值微积分与方程数值求解一、实验目的掌握求数值导数和数值积分的方法掌握代数方程数值求解的方法。掌握常微分方程数值求解的方法。二、实验内容要求：命令手工()输入1.求函数在指定点的数值导数。xx2f(x)二12x02x33x2,x=1,2,36x解：需要建立M脚本文件：clcclearx=1.0i=1f=inline('det([xxA2xA3;12*x3*xA2;026*x])')%内联函数whilex<=3.01为g增加一列做准备%x从1.0开始，以0.01的步长增加，一直到为g增加一列做准备%x从1.0开始，以0.01的步长增加，一直到3.01为止。x=x+0.01;endgt=1:0.01:3.01;%差分法近似求导%差分法近似求导%x=1的数值导数%x=2的数值导数%x=3的数值导数，即dx的第201个(最后一个)元素值。f1=dx(1)f2=dx(101)f3=dx(length(g)-1)运行结果：运行结果：2.用数值方法求定积分。⑴I=j2气cost2+4sin(21)2+1dt的近似值。10I=J2代ln(1+X)dt2o1+x2解：M文件：clc;clear;f=inline('sqrt(cos(t42)+4*sin(2*t)42+1)')%内联函数l1=quad(f,0,2*pi)%quad()函数是高精度的数值积分。g=inline('log(1+x)./(1+x42)')l2=quad(g,0,2*pi)3.分别用三种不同的数值方法解线性方程组。6x+5y-2z+5u=-49x-y+4z-u=13<3x+4y+2z-2u=13x-9y+2u=11解：M文件：clc;clear;A=[65-25;9-14-1;342-2;3-902]b=[-413111]'x=A\b%第一种方法，左除y=inv(A)*b%第二种方法，求逆矩阵A-诟，乘b[L,U]=lu(A)%第三种方法，矩阵的三角分解z=U\(L\b)4.求非齐次线性方程组的通解。2x+7x+3x+x=61234<3x+5x+2x+2x=412349x+4x+x+7x=21234解：先建立M函数文件，然后命令窗口中写命令。%%此部分需要建立M函数文function[x,y]=line_solution(A,b)[m,n]=size(A);y=[]；ifnorm(b)>0%非齐次方程组，norm()函数是计算矩阵范数ifrank(A)==rank([A,b])ifrank(A)==ndisp('有唯一解x');x=A\b;elsedisp('有无穷个解，特解x,基础解系y');x=A\b;y=null(Ay);%%null是用来求齐次线性方程组的基础解系的，如果加上'r',则求出的是一组最小正整数解。如果不加，则求出的是解空间的规范正交基。endelsedisp('无解')；x=[]；endelse%齐次方程组disp('有零解x');x=zeros(n,1);ifrank(A)vndisp('有无穷个解，基础解系y');y=null(A,'r');endend%%以下命令在窗口中输入clc;clear;formatrat%数据格式设置为有理数型A=[2731;3522;9417]b=[642]'[x,y]=line_solution(A,b)运彳丁结果：有无穷个解，特解x,基础解系yWarning:Rankdeficient,rank=2,tol=8.6112e-015.>Inline_solutionat11x=-2/1110/1100y=TOC\o"1-5"\h\z1/11-9/11-5/111/111001所以原方程组的通解是：-1/11-「-9/1「「-2/1「-5/111/1110/11X=k1+k0+0，其中k,k为任意常数12120105.求代数方程的数值解。(1)3x+sinx-ex=0在x0=1.5附近的根。(2)在给定的初值x0=1，y0=1，z0=1下，求方程组的数值解sinx+y2+Inz一7=0<3x+2y一z3+1=0x+y+z一5=0解：先建立M函数文件，然后命令窗口中写命令。方法一：clcclearf=inline('3*x+sin(x)-exp(x)')fzero(f,1.5)或者方法二：functiong=f(x)g=3*x+sin(x)-exp(x);命令窗口中输入：fzero('f',1.5)结果是：ans=1289/682(2).先建立M函数文件，然后命令窗口中写命令。%%此部分需要建立M函数文functionF=myfun(X)x=X(1);y=X(2);z=X(3);F(1)=sin(x)+yA2+log(z)-7;F(2)=3*x+2-zA3+1;F(3)=x+y+z-5;%%以下命令在窗口中输入X=fsolve('myfun',[1,1,1]',optimset('Display','off'))6.求函数在指定区间的极值。x3+cosx+xlogx(1)f(x)=在(0,1)内的最小值。ex(2)f(x,x)二2x3+4xx3一10xx+x2在[0,0]附近的最小值点和最小值。12112122解：(1Tclc;clear;formatlongf=inline('(xA3+cos(x)+x*log(x))/exp(x)');[x,fmin1]=fminbnd(f,0,1)(2)%先建立M函数文件％functionf=g(u)x=u(1);y=u(2);f=2*x.A3+4*x.*yA3-10*x.*y+y.A2;%命令窗口中写命令。[U,fmin2]=fminsearch('g',[0,0])U=1.0015701353166810.833488282765738fmin2=-3.324088491954234(以下选作题，是微分方程的数值解)7.求微分方程的数值解。x在[1.0e-9,20]'空-5強+y=0dx2dx<y(0)=0y'(0)=0解：M文件：functionxdot=sys(x,y)xdot=[y(2);(5*y(2)-y(1))/x];clc;clear;x0=1.0e-9;xf=20;[x,y]=ode45('sys',[x0,xf],[00]);[x,y]运行结果：xans=y'y0.0000000.5000001.0000001.5000002.0000002.5000003.0000003.5000004.0000004.5000005.0000005.5000006.000000

6.5000007.0000007.5000008.0000008.5000009.0000009.50000010.00000010.50000011.00000011.50000012.00000012.50000013.00000013.50000014.00000014.50000015.00000015.50000016.00000016.50000017.00000017.50000018.00000018.50000019.00000019.50000020.0000008.求微分方程组的数值解，并绘制解的曲线。〜y'=yyTOC\o"1-5"\h\z123y'=-yy〈213y'=—0.51yy312y(0)=0,y(0)=1,y(0)=1123解：令y1=x,y2=y,y3=z;这样方程变为:'x'=yz，自变量是ty'=—xz，自变量是tVz'=—0.51xy、x(0)=0,y(0)=1,z(0)=1M文件：functionxdot=sys(x,y)xdot=[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)];clc;clear;t0=0;tf=8;[x,y]=ode23('sys',[t0,tf],[0,1,1])plot(x,y)运行结果：x=0.00010.00050.00250.01250.06250.16320.30330.48290.71620.98491.26101.56781.95502.32872.70243.01533.29213.48893.64523.75383.86243.99414.16454.38354.65374.92655.22455.58616.03026.34286.6555

6.93717.15417.32387.45027.57657.70427.87068.000001.00001.00000.00011.00001.00000.00051.00001.00000.00251.00001.00000.01250.99991.00000.06240.99800.99900.16210.98680.99330.29650.95500.97730.45630.88980.94540.63500.77220.89120.79440.60690.82330.90690.42030.76170.97800.20660.71550.9975-0.06440.70160.9450-0.32580.73770.8127-0.58170.81410.6303-0.77550.89270.4130-0.90980.95520.2324-0.97160.98580.0795-0.99580.9980-0.0289-0.99860.9994-0.1367-0.98960.9949-0.2640-0.96340.9817-0.4187-0.90690.9538-0.5935-0.80340.9053-0.7644-0.64270.8373-0.8859-0.46090.7738-0.9656-0.25420.7235-0.9985-0.00140