例5.4

例5.4

#include<stdio.h>
#define SUM 100000
int main()
{
float amount,aver,total;
int i;
for (i=1;total=0;i<=1000;i++)
{
printf("please enter amount:");
scabf("%f",&amount);
total = total+amount;
if(total>=SUM) break;
}
aver=total/i;
printf("num=%d\naver=%10.2f\n",i,aver);
return 0;
}

在这里插入图片描述

### 实现四阶经典Runge-Kutta方法并比较 为了在MATLAB中实现四阶经典Runge-Kutta (RK4) 方法来求解微晶片固体激光器的速率方程组,并将其结果与 `ode43` 函数的结果进行比较,可以按照以下方式进行。 #### 定义速率方程组 假设微晶片固体激光器的速率方程组定义如下: ```matlab function dydt = laser_rate_eq(t, y) % 微晶片固体激光器的速率方程组 Ns = y(1); % 受激态粒子数密度 Pp = y(2); % 泵浦光功率 Pl = y(3); % 激光输出功率 gamma = 0.5; % 自发辐射率 tau_p = 1e-6; % 泵浦寿命 eta = 0.8; % 转换效率 alpha = 0.1; % 增益系数 dydt = zeros(3, 1); dydt(1) = -gamma * Ns + eta * Pp; dydt(2) = -(Pp / tau_p) + alpha * Ns * Pl; dydt(3) = alpha * Ns * Pp - Pl / tau_p; end ``` 此函数描述了受激态粒子数密度 \(N_s\)、泵浦光功率 \(P_p\) 和激光输出功率 \(P_l\) 的变化规律[^3]。 #### 编写四阶Runge-Kutta求解器 接下来,在MATLAB中编写一个通用的四阶Runge-Kutta求解器: ```matlab function [t, y] = runge_kutta_4(f, tspan, y0, h) % f 是导数函数句柄 % tspan 是时间区间 [t_start, t_end] % y0 是初值条件 % h 是步长 t_start = tspan(1); t_end = tspan(2); n = floor((t_end - t_start)/h); t = linspace(t_start, t_end, n+1)'; y = zeros(length(y0), length(t)); y(:,1) = y0'; for i = 1:n k1 = f(t(i), y(:,i)); k2 = f(t(i)+h/2, y(:,i)+(h*k1)/2); k3 = f(t(i)+h/2, y(:,i)+(h*k2)/2); k4 = f(t(i)+h, y(:,i)+h*k3); y(:,i+1) = y(:,i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4)'; end end ``` 这段代码实现了标准的四阶Runge-Kutta积分器,适用于任何给定的一阶常微分方程组[^1]。 #### 使用内置ODE solver (`ode43`) 对于对比实验部分,使用MATLAB自带的 `ode43` 来获得另一套解决方案作为参照。注意实际应用中应为 `ode45` 或者其他更合适的ODE求解器,这里遵照提问中的具体要求选用 `ode43` 进行说明: ```matlab [t_ode43, sol_ode43] = ode43(@laser_rate_eq, [0 15], [1; 1; 3]); subplot(1, 2, 1); plot(t_ode43, sol_ode43); title('ode43 解决方案'); xlabel('Time'); ylabel('Values'); legend({'Ns', 'Pp', 'Pl'}); ``` #### 结果可视化 最后一步是运行自定义的RK4算法并与上述得到的数据一起绘制图表来进行直观上的比较: ```matlab [h_t, h_y] = runge_kutta_4(@laser_rate_eq, [0 15], [1; 1; 3], 0.25); subplot(1, 2, 2); plot(h_t, h_y'); title('自编 Runge-Kutta 4 解决方案'); xlabel('Time'); ylabel('Values'); legend({'Ns', 'Pp', 'Pl'}); figure(); hold on; plot(t_ode43, sol_ode43, '--o'); plot(h_t, h_y', '-*r'); title('两种方法的对比'); xlabel('Time'); ylabel('Values'); legend({'ode43 Ns', 'ode43 Pp', 'ode43 Pl', ... 'RK4 Ns', 'RK4 Pp', 'RK4 Pl'},... 'Location','BestOutsidePlot'); grid on; hold off; ``` 通过这种方式可以在同一图上展示两个不同求解器产生的轨迹,从而便于观察两者之间的差异以及各自的特性。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值