DHU Matlab Experiment【5】作业记录_第五章、第六章

鱼犬

已于 2024-11-21 18:33:06 修改

阅读量1.1w

点赞数 18

分类专栏：学校课程记录文章标签： matlab 开发语言矩阵

于 2022-01-09 23:58:25 首次发布

本文链接：https://blog.csdn.net/weixin_45945586/article/details/122384408

版权

学校课程记录专栏收录该内容

28 篇文章

订阅专栏

这篇博客记录了使用Matlab进行数值微积分和解常微分方程的实验，涉及梯形法、Euler法、辛普生积分法等。内容涵盖河床断面积计算、椭圆周长积分、摩托车滑行距离、生产线停产时机、微分方程的数值解等实际问题的求解。

摘要由CSDN通过智能技术生成

写在前面

本博客用于记录（或者说是用来备份）我在2021a 高等数学实验课梁志勇老师的课上写的程序习题
所有的代码均通过编译，matlab版本为R2016a
课程链接：高等数学实验
由于题量的关系和对于篇幅的考虑，故一章到两章会做一篇博文，更多后续章节的答案可以点击我的头像查看或者点击下面的链接查看：
第一、第二章
 第三、第四章
 第七章
 网上参考答案

5.2.4 数值微积分

某河床的横断面如图5.8所示，为了计算最大的排洪量，需要计算它的断面积，试根据图示测量数据（单位：米）用梯形法计算其断面积。

clear;clc;
x = [0 4 10 12 15 22 28 34 40];
y = [0 1 3 6 8 9 5 3 0];
trapz(x,y)

在这里插入图片描述

求图5.8各测量点的坡度

clear;clc;
x = [0 4 10 12 15 22 28 34 40];
y = [0 1 3 6 8 9 5 3 0];
diff(y)./diff(x)

在这里插入图片描述

clear;clc;
t = 0:0.01:1.5;
x = log(cos(t));
y = cos(t)-t.*sin(t);
dydx=gradient(y,x)
[x,id] = min(abs(x-(-1)));%x最接近-1的编址
dydx(id)

在这里插入图片描述

求第（2）、（6）、（7）题积分的数值解

（2）

（6）注意：这里用了两种函数命名方式和其对应的求重积分函数

（7）

clear;clc;
fun = @(x,y) 1+x+y.^2;
ymin = @(x)-sqrt(2*x-x.^2);
ymax = @(x)sqrt(2*x-x.^2);
integral2(fun,0,2,ymin,ymax)

在这里插入图片描述

(椭园的周长) 用积分法计算下列椭园的周长

【法一】以dθ为底，做梯形积分

clear;clc;
th = 0:0.01:2*pi;
x = 2*cos(th);
y = 3*sin(th);
dx = gradient(x,th);
dy = gradient(y,th);
f = sqrt(dx.^2+dy.^2);
trapz(th,f)

在这里插入图片描述
【法二】用θ代换积分方程，做高精度数值积分
原理：

clc;clear;
f1 = @(th) 4*sqrt(4*(1-cos(th*pi/180).^2)+9*cos(th*pi/180).^2);
quadl(f1,0,pi/2)

在这里插入图片描述
【法三】用x代y，做高精度数值积分

clc;clear;
s=@(x)(1 - (9*x.^2)./(16*(x.^2/4 - 1))).^(1/2);
2*quad(s,-2,2)

在这里插入图片描述
由于高精度数值计算的原理未知，无法分析三种方法导致结果的差异性

5.3.3 数值微积分

(1) 用程序deriv.m求f(x)=x^2* sin(x^2 +3x-4)在x=1.3和x=1.5的导数，使精度达到10^-3。
(2) 编写用公式(5.20)求函数在某一点二阶导数达到指定精度的算法程序，并用此程序求f(x)=x^2* sin(x^2 -x-2)在x=1.4的二阶导数，使精度达到10^-3。

(1）将deriv.m当做函数使用即可

clear;clc;
fun = @(x) x.^2*sin(x.^2+3*x-4);
d1=deriv(fun,1.3,1e-6,1e-3)
d2=deriv(fun,1.5,1e-6,1e-3)

在这里插入图片描述
（2）
函数：

function d=diev(fname,a,h0,e)
h =h0;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/h^2;
d0 = d+2*e;
while abs(d-d0)>e
    d0 =d;h0=h;h=h0/2;
    d = (feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/h^2;
end
end

主程序：

clear;clc;
fun = @(x) x.^2.*sin(x.^2-x-2);
d=diev(fun,1.4,0.1,1e-3)

(辛普生积分法)编制一个定步长辛普生法数值积分程序。计算公式为

其中n为偶数，h=(b-a)/n,fi=f(a+(i-1)h). 并取n=50，应用于解习题5(1)。

n=2时的原理：（来源：陆老师的计算物理课）

函数：

function I=simp(fname,a,b,n)
%a:下限，b：上限
x = linspace(a,b,n+1);
y = feval(fname,x);
h = (b-a)/n;I = h/3*(y(1)+y(n+1)+2*sum(y(3:2:n))+4*sum(y(2:2:n)));
end

主程序及结果：
在这里插入图片描述

5.4.3 奶油蛋糕、作案时间

(摩托车)一个重5400kg的摩托车在以速度v=30m/s行驶时突然熄火，设滑行方程为

x为滑行距离，计算要滑行多长距离后, 速度可降至15m/s

移项积分即可，注意上下限，上限为15（末态），下限为30（初态）
在这里插入图片描述

(停产时间)某公司投资2000万元建成一条生产线。投产后，在时刻t 的追加成本和追加收益分别为

试确定该生产线在何时停产可获最大利润？最大利润是多少？

clear;clc;
fun = @(t) 17-t.^(2/3)-5-2.*t.^(2/3);
fplot(fun,[0,10]);%通过观察图像得知 函数单调递减，在8附近降为0
t = fzero(fun,8)%应在第八年停产
quadl(fun,0,t)-20%最大利润

在这里插入图片描述
即应在第8年停产，最大利润为1840万元

6.2.5 解常微分方程

解下列微分方程。
(1) y’=x+y, y(0)=1, 0<x<3(要求输出x=1, 2, 3点的y值)
(2) -x’=2x+3y, y’=2x+y,x(0)=-2.7,y(0)=2.8, 0<t<10, 作相平面图。
(3) y’'－0.01(y’)2+2y=sin(t), y(0)=0, y’(0)=1,0<t<5, 作y的图。

（1）
在这里插入图片描述
（2）

clear;clc;
fun=@(t,y)[-2*y(1)-3*y(2);2*y(1)+y(2)];
[t,y]=ode45(fun,[0 10],[-2.7;2.8]);
plot(y(:,1),y(:,2))

在这里插入图片描述
（3）
高阶导数y’‘化为一阶(y’)‘，多变量(y,y’)化为单变量x.
令x(1)=y,x(2)=y’,则有方程组x(1)‘=x(2),x(2)’=0.01x(2)^2-2x(1)+sin(t)
其初始值x(1)=0,x(2)=1.

clear;clc;
fun=@(t,x)[x(2);0.01*x(2)^2-2*x(1)+sin(t)];%fun表示两个方程的右端，第一个x(2)表示x(1)的导函数
[t,x]=ode45(fun,[0 5],[0;1]);x(end,1)
plot(t,x(:,1))

在这里插入图片描述

（解的“爆炸”）求一通过原点的曲线，它在(x,y)处的切线斜率等于2x+y^2,0。若x上界增为1.58,1.60会发生什么？

%由题意，易知y'=2x+y^2,y(0)=0
clear;
fun=inline('2*x+y^2','x','y'); 
[x,y]=ode45(fun,[0 1.57],0)
%[x,y]=ode45(fun,[0 1.58],0)
%[x,y]=ode45(fun,[0 1.60],0)
plot(x,y)

若上界增加，函数爆炸，t步长要求过小
在这里插入图片描述

6.3.2 Euler法和刚性方程组

求解刚性方程组

函数

function f=eg6_7f(x,y)
f(1)=-1000.25*y(1)+999.75*y(2)+0.5;
f(2)=999.75*y(1)-1000.25*y(2)+0.5;
f=f(:)
end

主程序

clear;clc;
[x,y]=ode45(@eg6_7f,[0 5],[1;-1]);

在这里插入图片描述

6.4.4 导弹系统的改进

(温度过程)夏天把开有空调的室内一支读数为20℃的温度计放到户外，10分钟后读25.2℃, 再过10分钟后读数28.32℃。建立一个较合理的模型来推算户外温度。(提示：用牛顿冷却定律(Newton’s law of cooling)）

%由牛顿冷却定律：T'=k(c-T),T(0)=20
clear;clc;
dsolve('DT=a*(c-T)','T(0)=20','t');%Ch7
%用上式得到解析解c - exp(-a*t)*(c - 20)，下面利用T(10)=25.2, T(20)=28.32拟合
fun=@(c,t)c(1) - exp(-c(2)*t)*(c(1) - 20);
lsqcurvefit(fun,[30 1],[10 20],[25.2 28.32])%初始值c=30是估计的，k=1是一个尝试值(根据建模意义知k>0)

在这里插入图片描述
解得户外温度c=33,比例系数k=0.05.

(广告效应)某公司生产一种耐用消费品，市场占有率为5%时开始做广告，一段时间的市场跟踪调查后，该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比，且估得此比例系数为0.5。
(1) 建立该问题的数学模型，分别求其解析解和数值解，并作比较；
(2) 厂家问：要做多少时间广告，可使市场购买率达到80%？

（1）

%根据题设，列出方程：x'/x=0.5*(1-x)  x'=0.5*x*(1-x) 
clear;clc;
dsolve('Dx=0.5*x*(1-x)','x(0)=0.05','t');%解析解：1/(exp(log(19) - t/2) + 1)
fun=@(t,x)0.5*(1-x)*x; [t,x]=ode45(fun,[0 10],0.05);%数值解
%作图比较
fplot(@(t)1./(1+19*exp(-1/2.*t)),[0,10]);hold on;
plot(t,x,'o');hold off; 
%易知数值解与解析解基本重合

在这里插入图片描述
（2）

clear;clc;
a=@(t)1./(1+19*exp(-1/2.*t))-0.8;%t:所需时间
fzero(a,[0 10])

在这里插入图片描述