数学实验课MATLAB实验报告二（题目+代码）

前言

2022年10月21日晴转多云转晴然后黑天了，不冷。今天有一件要紧的事要做，但我就是要先写完这个再去做。

1.(1)

题目

解微分方程
$$\{\begin{array}{l}\frac{d^{3}y}{d{x}^3}=(\frac{d^{2}y}{d{x}^2}-1{)}^2-\frac{dy}{dx}-y^2,\\ y(0)=0,y^{\prime }(0)=1,y^{\prime \prime }(0)=-1.\end{array}$$
（考虑$x\in [0,20]$，画出数值解及其一阶、二阶导数曲线图形）.

代码

高阶微分方程，先化成一阶微分方程组，再求解。
令$y_1=y,y_2=y^{\prime },y_3=y^{\prime \prime }$，则原方程化为方程组
$$\{\begin{array}{l}{y}_1^{\prime }=y_2,\\ y_2^{\prime }=y_3,\\ y_3^{\prime }=(y_3-1{)}^2-y_2-y_1^2,\\ y_1(0)=0,y_2(0)=1,y_3(0)=-1.\end{array}$$
所以在编写函数脚本文件时，要写成下面这样。

function dy=f1(t,y)
dy=[y(2);y(3);(y(3)-1).^2-y(2)-y(1).^2];
end

将上面的代码保存为f1.m。其中y是一个三列的向量，分量1为$y_1$，分量2为$y_2$，分量3为$y_3$。下面我们就要在ode45中调用f1啦：

[T,y]=ode45('f1',[0 20],[0 1 -1]);
plot(T,y(:,1),'-',T,y(:,2),'o',T,y(:,3),'x')

代码中[0 20]中的0是初值条件中的初值点，后面的[0 1 -1]是初值条件。画出的图像如下：

1.(2)

题目

求
$$\{\begin{array}{l}\frac{dy}{dx}=-0.01y-99.99z,\\ \frac{dz}{dx}=-100z,\\ y(0)=2,z(0)=1.\end{array}$$
的数值解，并画出数值解的图像。

代码

一样的，令$y_1=y,y_2=z$，则有
$$\{\begin{array}{l}{y}_1^{\prime }=-0.01y_1-99.99y_2,\\ y_2^{\prime }=-100y_2,\\ y_1(0)=2,y_2(0)=1.\end{array}$$

function dy=f2(t,y)
dy=[-0.01*y(1)-99.99*y(2);-100*y(2)];
end

将上面这个代码保存为f2.m。调用：

[T,y]=ode45('f2',[0 500],[2 1]);
plot(T,y(:,1),'-',T,y(:,2),'.')

刚开始画的区间是0到20，但是很丑，很怪，不断调整区间，最终觉得0到500画出来的图还能入眼：

1.(3)

题目

求下列方程的通解及特解。
$$\{\begin{array}{l}{x}^2{y}^{\prime \prime }+x{y}^{\prime }+(x^2-\frac{1}{4})y=0,\\ y(\frac{\pi }{2})=2,y^{\prime }(\frac{\pi }{2})=-\frac{2}{\pi }.\end{array}$$

代码

%求通解
y=dsolve('x.^2*D2y+x*Dy+(x.^2-1/4)*y=0','x');
disp('通解')
disp('y=')
pretty(y)
%求特解
y0=dsolve('x.^2*D2y+x*Dy+(x.^2-1/4)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x');
disp('特解')
disp('y0=')
pretty(y0)

用dsolve求解。求通解不需要写初值，求特解时要写初值。结果如下：

1.(4)

题目

已知$y(0)=-0.2,y^{\prime }(0)=-0.7$用数值方法求$VanderPol$方程$y^{\prime \prime }+\mu (y^2-1)y^{\prime }+y=0$的解，并画出当$\mu =1000$时的图像。

代码

倦了。。。令$y_1=y,y_2=y^{\prime }$。So easy!!!

function dy=vdp1000(t,y)
dy=[y(2);-1000*(y(1).^2-1)*y(2)-y(1)];
end

保存为vdp1000.m。。。

[T,Y]=ode45('vdp1000',[0 0.05],[-0.2;-0.7]);
plot(T,Y(:,1),'-');
title('$Solution\ of\ van\ der\ Pol\ Equation,\mu=1000$','interpreter','latex');

哦，这里注意一下，最开始我给的区间是0到20，画出来的图巨丑，后来发现它就这么点好看的地方，我就把区间改成0到0.05咯。结果：

2.

题目

凶杀案作案时间问题：受害者的尸体于晚上7:30被发现，法医于晚上8:20赶到凶案现场，测得尸体温度为32.6℃；一小时后，当尸体即将被抬走时，测得尸体温度为31.4℃，室温在几个小时内始终保持21.1℃。此案最大的嫌疑犯张某声称自己是无罪的，并有证人说：“下午张某一直在办公室上班，5:00时打完电话后就离开了办公室”。从张某到受害者家（凶案现场）步行需5分钟，现在的问题是，张某不在凶案现场的证言能否被采信，使他排除在嫌疑犯之外。（提示：1、Newton冷却定理告诉我们“物体在介质中冷却速度同该物体温度与介质温度之差成正比” 2、自己假设人体的正常体温的温度值；3、求解非线性方程的MATLAB命令有fzero,solve）

代码

这个题目我们要求的是受害者死亡的时刻 $t_0$。
首先介绍一下，%%在脚本文件中是分节符。在第一节，我们把体温随时间变化的关系式解出来了。将第一次测温的8:20作为0时刻，求解即可。运行第一节，我们得到关系式为 $T=\frac{23e^{kt}}{2}+\frac{211}{10}$
在第二节中，使用一小时后9:20的条件解出 $k$ 的值。这里选择时间 $t$ 的单位为小时，故 $t=1$，用solve解非线性方程即可，解出 $k$ 的值赋给 $k_0$ 。
下面利用得到的特解 $T=\frac{23e^{k_{0}t}}{2}+\frac{211}{10}$求温度 $T=37$ 的时间 $t_0$。为什么体温是37呢？因为人受到惊吓时体温会上升那么一丢丢。

%%求解温度变化曲线
clc,clear
T=dsolve('DT-k*(T-21.1)=0','T(0)=32.6','t')
%%
syms k
k0=solve((23*exp(k))/2 + 211/10-31.4==0);
syms t
t0=solve((23*exp(k0*t))/2 + 211/10-37==0);
time=8+1/3-ceil(abs(t0))+mod(t0,1);
a=sprintf('%d',floor(time));%时
b=sprintf('%d',floor(mod(time,1)*60));%分钟
disp(['案发时间为' a ':' b '左右'])

代码到time=…之前就已经把题目求解完了，后面只是将输出结果美化。blingbling美美哒的结果：

心机之蛙一直摸你肚子！案发时间是5:23左右，那么张某的嫌疑就无法排除了。

总结

妹想到这么快又摸MATLAB了…心碎。随便做了做，若有错误请指正（抱拳）（充满侠气地离开）。