My Numerical Solution doesn't align with the exact solution

4 次查看（过去 30 天）
Prince Nino
Prince Nino 2024-9-19
评论： Prince Nino 2024-9-21
Ran in:
This is a 2nd order differential equation and I'm using ODE45 for my numerical Solution. My Numerical Solution is way too high. Notice that the values for numerical is *10^8. Additionally, I first solve my constants in order to make my equation code at matlab simpler.
My initial conditions are vc(0) = 11; vc' = 77459.66692
Exact Solution:
clear all
close all
clc
C = 3*10^(-6);
R = 12;
L = 2*10^(-3);
a = 0;
B = (6/(sqrt( L*C )))/(sqrt((1/(L*C))-(R/(2*L))^2));
Vc_0 = [6/(sqrt( L*C )) , 11];
odefun = @(t,Vc) [Vc(2);...
((-600)*(Vc(2)))-((50000000/3)*(Vc(1)))];
tspan = linspace(0, 0.04, 4000000);
[t,Vc] = ode45(odefun, tspan, Vc_0);
subplot(2,1,1);
plot(t,Vc)
xlabel('time(s)', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
ylabel('Vc', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
title('Numerical Methods', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
texact=[0:0.00000001:0.04];
Vexact=exp((-300)*texact).*(a*cos(((sqrt(596760000)*texact)/6))+B*sin(((sqrt(596760000)*texact)/6)));
subplot(2,1,2);
plot(texact,Vexact)
xlabel('time(s)', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
ylabel('VCexact', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
title('Exact Solution', 'FontSize',16, 'FontName','Arial','FontWeight','bold')
  3 个评论
显示 1更早的评论
Prince Nino
Prince Nino 2024-9-20
That's what I was wondering. Why would our instructor have initial values for our exact solution be different from the numerical solution.
I have tried inputting the C,R,L Variables but it seems it gives me different answers(idk why). I've tried solving it manually. So I opted to solve it when it is already numerically valued rather than being it variable.
And yes, all the units are correct.
Sam Chak
Sam Chak 2024-9-20
Ran in:
@Prince Nino, Because you didn't use the values derived from the parameters R, L, C in the code. Instead, you entered the coefficients '600' and '50000000/3' manually.
%% Parameters
C = 3*10^(-6);
R = 12;
L = 2*10^(-3);
R/L
ans = 6000
1/(L*C)
ans = 1.6667e+08
50000000/3
ans = 1.6667e+07

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

James Tursa
James Tursa 2024-9-19
编辑：James Tursa 2024-9-19
Assuming your state is [vc,vc'] in that order, then it appears your initial conditions are backwards and should be this instead:
Vc_0 = [11, 6/(sqrt( L*C ))];
I haven't looked any deeper into your problem yet ...
Sam Chak
Sam Chak 2024-9-20
Ran in:
Hi @Prince Nino
Perhaps you should verify the results with your Professor. You can also change the initial condition.
%% Numerical Solution
C = 3*10^(-6);
R = 12;
L = 2*10^(-3);
f = @(t, x) [ x(2);
-(R/L)*x(2) - 1/(L*C)*x(1)];
tt = linspace(0, 0.004, 4001);
ic = [0; 6/sqrt(L*C)]; % initial condition
[t, x] = ode45(f, tt, ic);
%% Analytical Solution
a = 1;
b = R/L;
c = 1/(L*C);
rpart = b/(2*a);
ipart = sqrt(c/a - rpart^2);
alpha = ic(1);
beta = (ic(1) + ic(2))/ipart;
xt = exp(- rpart*tt).*(alpha*cos(ipart*tt) + beta*sin(ipart*tt));
%% Plot results
tl = tiledlayout(2, 1);
nexttile
plot(t, x(:,1)), grid on, title('Numerical Solution')
nexttile
plot(tt, xt), grid on, title('Analytical Solution')
xlabel(tl, 'Time')
ylabel(tl, 'x(t)')

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by