Simpson 法求数值微分

for i=2:181 th1=theta1(i);

theta3(i)=fzero(‘ex714f’,theta3(i-1));%调用fzero函数求theta3 end

subplot(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 %画角速度曲线 3.子程序（函数程序）ex714f.m function y=ex714f(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;

在上述程序中主要th1是一个标量，而theta1在ex714b中是一个数组，因为函数ex714f中用到的是特定点的角度，是一个标量，所以不能用theta1,引入了th1作为其当前量。 程序运行结果

输入输出角位置关系2.62.42.223ateht1.81.61.402468theta1输出角速度6040200-20-40-60-8002468theta1w3

实验7用电压V=10伏的电池给电容器充电，电容器上t时刻的电压V（t）=V-（V-V0）exp(-t/τ),其中V0是电容器的初始电压，τ是充电常数，试由下面一组t，V数据确定V0和τ。 t/s 0.5 V/V 6.36 分析： 1 6.48 2 7.26 3 8.22 4 8.66 5 8.99 7 9.43 9 9.63 采用最小二乘法进行拟合，对V（t）=V-（V-V0）exp（-t/τ）两边求自然对数得到：log(V-V(t))=log(V- V0)-t/τ

令k=-1/τ，w=log（V- V0），x=t，y=log（V-V（t））。 得到方程：y=kx+w 流程图：

结束 输出τ、V0以及相应的函数图象 得到τ=1/p(1), V0=10-exp(p(2)) 利用多项式函数进行拟合：[p, s]=polyfit(t, V1, 1)进行一次多项式拟合。 对原函数V（t）=10 -（10-V0）exp （-t/τ）两边取自然对数，可得：log(10-V(t))=log(10- V0)-t/τ，变形后可得：w=log（V- V0），x=t，y=log（V-V（t）），k=-1/τ 开始 程序代码： clear;clc;

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]; V1=log(10-V); [p, s]=polyfit(t, V1, 1); [V1p, delta]=polyval(p, t, s); t0= -1/p(1) V0 = 10 - exp(p(2))

plot(t, V, '--r', t, 10 - (10 - V0)*exp(-t/ t0)) legend('实验数据','拟合数据') 结果： t0 =

3.5269e+000 V0 =

5.6221e+000