画orr-sommerefeld 问题的中性曲线

*若价格不公道，可以让提问者在平台追加赏金哦，平台是您利益的保证

已完成

134****0008 2年前发布

悬赏：300.0 元

问题详情：分类: 其他-Matlab

1. orr-sommerefeld 问题是一个常微分方程特征值问题。我会给一个代码用于求特征值c(Matlab的Chebfun包运行)，在提供的代码中我任意给定了R和ak（也就是）的数值，这组R和对应很多个特征值c，我们只需要考虑最大的那个c（代码中为first_eigenvalues_c)。在提供的代码中，我给定R=10000，=1，则输出的first_eigenvalue_c=0.237528104235363 + 0.003737973691874i，我们可以看到它是一个复数。任务：找到多组R和，使得first_eigenvalue_c的虚部=0（我们无法让它真的=0，让它的虚部<0.000001），连接这些满足条件的R和得到中性曲线（以R和为坐标）。我已经用笨方法得到大致曲线如下图，但是运行时间太长，希望大神可以用牛顿迭代或者其他方法缩短运行时间

orr-sommerefeld equation:
(U-c)(v''-^2*v)-vU''=-i/(*R)*(v''''-2^2v''+^4v)
U=1-y^2
v=0 at y=1,y=-1
v'=0 at y=1,y=-1

clear all; clc;
Re = 10000;
ak = 1;

B = chebop(-1,1);
B.op = @(x,u) diff(u,2) - ak^2*u;
A = chebop(-1,1);
A.op = @(x,u) (diff(u,4)-2*ak^2*diff(u,2)+ak^4*u)/Re - ...
    2i*ak*u - 1i*ak*(1-x^2)*(diff(u,2)-ak^2*u);
A.lbc = @(u) [u; diff(u)];
A.rbc = @(u) [u; diff(u)];

[V,lam] = eigs(A,B,100,'LR');

%reorder eigenvalues
d=diag(lam);
[d1,iii] = sort(d,'descend','ComparisonMethod','real');
c = d1/ak*1i;
%eigenvalues we wanted

format long
first_eigenvalue_c = d1(1)/ak*1i