converting python code to matlab and getting it plotted
13 ビュー (過去 30 日間)
古いコメントを表示
I been trying to get this python code to work in matlab. Besides converting issues I probably have a lot of other mistakes in it as I am just a beginner.
% constants
L = 100; % length of bar in km
dx = 1; % spatial grid spacing km
beta = 4; % shear wave velocity in km/s
dt = 0.1; % time spacing in seconds
T = 5; % total time in sec
N = int(T/dt); % number of time steps
tmax = 30; % wave run time sec
t = [0:dt:tmax]; %time vec
nt = length (t);
% Define initial conditions
u = zeros(int(L/dx)+1); % displacement array
u(int(50e3/dx)) = 1; % source at u(50 km)
% intialize 3 displacement vectors
%u3 = new solution
%u2 = previous solution
%u1 = previous previous solution
% displacement arrays
u1 = zeroes (1, 101);
u2 = u1;
u3 = u1;
% boundary conditions
u(0) = 0;
u(-1) = 0;
% Define finite difference coefficients
A = (beta*dt/dx)^2
B = 2 - 2*A
% Solve wave equation using finite differences
for n = range(1,N)
u(1:-1) = B*u(1:-1) - u(1:-1) + A*(u (2: + u :-2))
u(int(50e3/dx)) = np.sin(2*np.pi*n*dt/5)
end
% Plot solution
figure(1); hold on; grid on; axis equal
xline(0,'-','linewidth',2); yline(0,'-','linewidth',2) %adding x and y line
set(0,'DefaultLineLineWidth',5,'DefaultAxesFontSize',12)
x = linspace(0, L, int(L/dx)+1)
plot(x, u)
xlabel('Distance (m)')
ylabel('Displacement')
show()
% initial condition
u3(50) = sin(pi*tlen/5)^2
% boundary conditions
u3(0) = u3(1) % stress-free boundary
u3(L) = 0 % fixed boundary
% Define plotting function
figure plot_u(u, t):
x = arange(0, L+dx, dx)
plot(x, u)
title('t = {t:.2f} s')
xlabel('x (km)')
ylabel('u (m)')
ylim(-0.1, 0.1)
show()
%%% Solve PDE using finite differences
t = 0
if t <= tlen
t = dt
for i [1, L]
rhs = beta^2 * (u2[i+1] - 2*u2[i] + u2[i-1]) / dx^2
u3[i] = 2*u2[i] - u1[i] + dt^2 * rhs
% Apply boundary conditions
u3(0) = u3(1) % stress-free boundary
u3(L) = 0 % fixed boundary
if t <= tlen
u3(50) = sin(np.pi*t/tlen)^2
% Update displacement arrays
u1 = copy(u2)
u2 = copy(u3)
% Plot every 4 s
if t % 4 < dt:
plot_u(u2, t)
end
end
end
end
%% Calculate velocities of pulses
t1 = 1.5 % time to measure velocity of first pulse
t2 = 12.5 % time to measure velocity of second pulse
x1 = 25 % distance to measure velocity of first pulse
x2 = 75 % distance to measure velocity of second pulse
vel1 = (u2(x1+1,t1) - u2(x1-1,t1)) / (2*dx*dt)
vel2 = (u2(x2+1,t2) - u2(x2-1,t2)) / (2*dx*dt)
fprintf('Velocity of first pulse: {vel1:.2f} km/s')
fprintf('Velocity of second pulse: {vel2:.2f} km/s')
figure(2) % displacement at endpoints as a function of time
x_endpoints = [0, L]
u_endpoints = u2*[0,':'], u2*[L,':']
for i % in range(2):
plot(arange(0, tlen+dt, dt), u_endpoints[i])
xlabel('t (s)')
ylabel('f u(x_endpoints[i] km) (m)')
show()
% Plot displacement at crossing
(initialize t, dx, dt, tlen, beta and u1, u2, u3 arrays)
while t <= 33
t = t + dt;
for i = 2:99
rhs = beta^2 * (u2(i+1) - 2*u2(i) + u2(i-1)) / dx^2;
u3(i) = dt^2 * rhs + 2*u2(i) - u1(i);
end
u3(1) = u3(2); % stress-free boundary
u3(100) = 0; % fixed boundary
if t <= tlen
u3(50) = sin(pi * t/tlen)^2;
end
u1 = u2;
u2 = u3;
% output u2 at desired times
end
3 件のコメント
採用された回答
Sarthak
2024 年 2 月 26 日
Hi Kathleen,
As per my understanding, there is no direct way to convert python code to MATLAB.
However you can try the following approaches:
- You have to directly write the code from scratch. You should be able to find direct MATLAB coutnerparts to PYTHON functions.
- You can use the MATLAB API to call python script in MATLAB: https://www.mathworks.com/help/matlab/call-python-libraries.html
- As of MATLAB 8.4(R2014b) you can call Python directly from MATLAB. For more information, please go through the following documentation: https://www.mathworks.com/help/matlab/python-language.html
- Can also leverage the following library on FileExhange: https://www.mathworks.com/matlabcentral/fileexchange/114455-python-for-matlab-development
Please refer to the following MATLAB answers post with similar issue: https://www.mathworks.com/matlabcentral/answers/416129-how-to-convert-python-code-into-matlab
I hope the above solutions successfully resolve your query!
0 件のコメント
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Call Python from MATLAB についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!