My Numerical Solution doesn't align with the exact solution
123 views (last 30 days)
Show older comments
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 Comments
Sam Chak
on 20 Sep 2024 at 12:09
@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
1/(L*C)
50000000/3
Answers (2)
James Tursa
on 19 Sep 2024 at 19:35
Edited: James Tursa
on 19 Sep 2024 at 19:38
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
on 20 Sep 2024 at 5:19
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)')
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!