My Numerical Solution doesn't align with the exact solution

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')

[Torsten]

What strikes me is that your exact solution does not start at (0,11). The fact that you choose alpha = 0 in the "analytical solution" automatically leads to Vc(0) = 0 which contradicts your initial value for Vc.

And the order of your initial conditions in the numerical solution is wrong:

Vc_0 = [6/(sqrt( L*C )) , 11];

I suggest that you insert the expressions depending on C, R and L for the coefficients of the ODE instead of error prone numerical values. E.g. exp((-300)*texact) must be exp(-3000*texact) instead (I didn't check the other numerical values).

Are you sure that the initial conditions for Vc and dVc/dt are given in SI units ?

[Prince Nino]

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]

@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

Answer 1

[James Tursa]

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 ...

[Prince Nino]

Isn't it [vc', vc] ?

Answer 2

[Sam Chak]

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)')

[Prince Nino]

Thanks