高斯扩散模型实现

35 篇文章 1 订阅
34 篇文章 8 订阅

算法公式来源于(高斯扩散模型-高斯烟羽大气污染扩散模型 - 简书 (jianshu.com)),感谢大哥的奉献,接下来我们在有算法由公式的情况下来实现。

一顿代码巴啦啦。。。。

1.原始数据(原始数据就是一个区域加一个污染源)

2.代码参数界面

3.结果

 

4.结论,值得注意的是,根据公式计算出的结果在靠近源点位置有几个异常大的值,需要剔除异常值。 

  • 2
    点赞
  • 12
    收藏
    觉得还不错? 一键收藏
  • 13
    评论
这里提供一个基于遗传算法的高斯烟雨扩散模型的 Matlab 代码,供参考。 ``` % 遗传算法参数 popSize = 100; % 种群大小 maxGen = 100; % 最大迭代次数 pc = 0.8; % 交叉概率 pm = 0.01; % 变异概率 % 模型参数 xRange = [0, 1000]; % x 范围 yRange = [0, 1000]; % y 范围 sigma = 100; % 高斯函数标准差 Q = 10; % 源强度 u = 2; % 风速 v = 0; % 风向 D = 0.1; % 扩散系数 dt = 1; % 时间步长 % 初始化种群 pop = rand(popSize, 2) .* [diff(xRange), diff(yRange)] + [xRange(1), yRange(1)]; % 迭代 for gen = 1:maxGen % 计算适应度 fitness = zeros(popSize, 1); for i = 1:popSize x = linspace(xRange(1), xRange(2), 100); y = linspace(yRange(1), yRange(2), 100); [X, Y] = meshgrid(x, y); z = Q / (2 * pi * u * sigma) * exp(-((X - pop(i, 1)).^2 + (Y - pop(i, 2)).^2) / (2 * sigma^2)) .* exp(-D*dt*(u*(X - pop(i, 1)) + v*(Y - pop(i, 2)))^2); fitness(i) = sum(sum(z)); end % 选择 [sortedFitness, indices] = sort(fitness, 'descend'); elite = pop(indices(1:ceil(popSize/10)), :); roulette = cumsum(sortedFitness) / sum(sortedFitness); parents = zeros(popSize, 2); for i = 1:popSize r = rand; for j = 1:popSize if r < roulette(j) parents(i, :) = pop(indices(j), :); break; end end end % 交叉 offspring = zeros(popSize, 2); for i = 1:popSize/2 if rand < pc idx1 = randi(popSize); idx2 = randi(popSize); alpha = rand; offspring(2*i-1, :) = alpha*parents(idx1, :) + (1-alpha)*parents(idx2, :); offspring(2*i, :) = alpha*parents(idx2, :) + (1-alpha)*parents(idx1, :); else offspring(2*i-1, :) = parents(randi(popSize), :); offspring(2*i, :) = parents(randi(popSize), :); end end % 变异 for i = 1:popSize if rand < pm offspring(i, 1) = offspring(i, 1) + randn * diff(xRange) / 10; offspring(i, 2) = offspring(i, 2) + randn * diff(yRange) / 10; end end % 更新种群 pop = [elite; offspring]; end % 计算最优解 x = linspace(xRange(1), xRange(2), 100); y = linspace(yRange(1), yRange(2), 100); [X, Y] = meshgrid(x, y); z = zeros(size(X)); for i = 1:popSize z = z + Q / (2 * pi * u * sigma) * exp(-((X - pop(i, 1)).^2 + (Y - pop(i, 2)).^2) / (2 * sigma^2)) .* exp(-D*dt*(u*(X - pop(i, 1)) + v*(Y - pop(i, 2)))^2); end z = z / popSize; % 绘图 surf(X, Y, z); xlabel('x'); ylabel('y'); zlabel('concentration'); ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 13
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值