材料科学与工程学院。
2024年4月12日。
一、编程实现“simpson法求数值微分”，并举一例应用之。
实例： 分别采用5，10，100个离散点的辛普森数值微分法求函数f=ex在x=2.5的导数值。
1.算法说明：
辛普森数值微分是用来求等距节点在节点处的导数的，辛普森数值微分公式如下：
其中，yn=f(xn),xn=x0+nh。
如果端点导数值-f’(x0)和-f’(xn)未知，则将它们用中点微分公式近似，这时的辛普森微分公式为：
辛普森数值法函数为：cisimpson；
功能：辛普森数值法求已知函数在某点的导数值；
调用格式：df=cisimpson(func,x0,n,h)。
2.流程图。
simpson数值微分法流程图：
实例求解流程图：
3.源程序**。
simpson数值微分法源程序：
function df=cisimpson(func,x0,n,h)
辛普森数值法求已知函数func在x0点的导数值。
函数名：func
求导点：x0
将已知函数离散的数据点数：n
离散步长：h
导数值：df
if nargin ==2以下是参数的判断过程。
h=0.1;
n=5;else
if (nargin ==3)
if(n<5)
disp('n不能小于5!')
return;
elseh=0.1;
endelse (nargin ==4 &&h ==0.0)
disp('h不能为0!')
return;
endend
for(i=1:n) %这个循环计算节点的函数值。
if (mod(n,2) =0)
y(i)=subs(sym(func),findsym(sym(func)),x0+(i-n/2)*h);
elsey(i)=subs(sym(func),findsym(sym(func)),x0+(i-(n+1)/2)*h);
endend
f(1)=(y(3)-y(1))/2*h);
f(2)=(y(n)-y(n-2))/2*h); 这两行用中心微分法给出端点的导数。
b(1:n-2,1) =zeros(n-2,1);
b(1,1)=3*(y(3)-y(1))/h-f(2);
b(n-2,1)=3*(y(n)-y(n-2))/h-f(2);
for(i=2:(n-3))
b(i,1) =3*(y(i+2)-y(i))/h;
end这一块是辛普森公式的右边的列向量。
for(i=1:n-2)
for(j=1:n-2)
if((i ==j+1) |j ==i+1))
a(i,j)=1;
else if(i ==j)
a(i,j) =4;
endend
endend这一块是系数矩阵。
[q,r]=qr(a);
df = r\(q\b用qr分解法求解。
if(mod(n,2) =0)
df = df(n/2);
elsedf = df((n+1)/2);
end这里是求出x0处的导数值。
实例求解源程序：
df = cisimpson('exp(x)',2.5) %采用默认的5个离散点。
df = cisimpson('exp(x)',2.5,10) %采用10个离散点计算。
df = cisimpson('exp(x)',2.5,100) %采用100个离散点计算。
4.程序运行结果。df =
df =df =
一、编程解决以下科学计算和工程实际问题。
1）[例7-1-4]四连杆机构如图7-1-4所示，输入杆l1的转角……并求其角速度和角加速度。（图略）
1.流程图。
2.源程序**。
global l0 l1 l2 l3 th1
l0=20;l1=8;l2=25;l3=20;
%输入基线及三根杆的长度l1、l2、l3
theta1=input('当前角theta1=')
theta3=input('对应于theta1的theta3近似值=')
th1=theta1;theta3=fzero('exn714f',theta3);
%求当前输出角theta3
theta2=asin((l3*sin(theta3)-l1*sin(theta1))/l2);
w1=input('w1=')
w3=l1*w1*cos(pi/2-theta1+theta2)/(l3*cos(theta3-pi/2-theta2));
global l0 l1 l2 l3 th1
l0=20;l1=8;l2=25;l3=20;
%输入基线及三根杆的长度l1、l2、l3
w1=input('杆1角速度w1=')
theta1=linspace(0,2*pi,181);
%把杆1每圈分为180份，间隔2度。
theta3=input('对应于theta1最小处的theta3(近似估计值)='
dt=2*pi/180/w1; %杆1转2度对应的时间增量。
th1=theta1(1);theta3(1)=fzero('exn714f',theta3);
%求初始输出theta3
for i=2:181
th1=theta1(i);
theta3(i)=fzero('exn714f',theta3(i-1));
%调用fzero函数逐次求theta
endsubplot(1,2,1),plot(theta1,theta3),ylabel('theta3'),grid;
%画曲线。w3=diff(theta3)/dt;
%求杆3的角速度，注意求导数后数组长度小于1
subplot(1,2,2),plot(theta1(2:length(theta1)),w3);grid
%画角速度曲线。
function y=exn714f(x)
global l0 l1 l2 l3 th1
y=l1.*cos(th1)+l2*sqrt(1-(l3*sin(x)-l1*sin(th1)).2/l2/l2)-l3*cos(x)-l0;
3.程序运行结果。
运行根据提示，输入w1=100，在theta1=0处，设theta3的近似初值为1弧度，得到曲线如图a所示，相应的角速度变化规律如图b所示。
2）实验7 用电压v=10伏的电池给电容器充电，电容器上时刻的电压……，是充电常数，试由下面一组t，v数据确定v0和τ。
1. 流程图。
2.源程序**。
function f=curvefun1(x,t)
f=10-(10-x(1))*exp(-t/x(2));
t=[0.5 1 2 3 4 5 7 9];
v=[6.36 6.48 7.26 8.22 8.66 8.99 9.43 9.63];
x0=[0.2 0.05];
x=lsqcurvefit('curvefun1',x0,t,v);
f=curvefun1(x,t);
plot(t,v,'*t,f,'g')
x3.程序运行结果。
运行得到：x=
MATLAB课程设计
1 求被控对象传递函数g s 的matlab描述。num 789 6312 11835 den 1 14 56 64 0 0 gs tf num,den transfer function 789 s 2 6312 s 11835 s 5 14 s 4 56 s 3 64 s 2 2 求被控对象脉冲...
MATLAB课程设计
课程设计。题目 matlab计算器。姓名 班级 学院 专业 完成时间。1总体设计。该计算器程序主要是matlab来制作，界面主要由四个静态文本框 21个运算按钮和两个动态文本框组成。实现的运算功能有四则运算 加 减 乘 除。而且添加了括号使人们使用时更加简单。这些计算功能主要调用了matlab的自定...
MATLAB课程设计
matlab课程设计。如图所示，为测量系统的示意图，它由两个能相互转动的连杆，角度编码器和滚轮等组成。o1为固定点，o2点为转动点，o3点为滚轮的中心，连杆的有效长度分别为l1和l2。任一位置时，连杆1相对于某基准位置的角度为 1,两连杆的相对角度为 2。其中对于 1，取垂直方向为基准线，在基准线左...